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ABSTRACT 


Most recent studies involving multiple swingby interplanetary tra- 
jectories have been made using a simplified model consisting of a se- 
quence of heliocentric conic arcs matched in relative hyperbolic exess 
velocity at each planetary encounter. This model provides adequate re- 
sults for preliminary mission planning and analysis but as more ad- 
vanced investigations are undertaken, an accurate N-body reference tra- 
jectory becomes necessary. This thesis presents a technique for the 
rapid determination of such a reference trajectory. 

The gap between the simple conic model and the integrated N-body 
trajectory is bridged in two steps. The first of these utilizes a 
model of the trajectory consisting of alternating planetocentric and 
heliocentric conic legs corresponding to trajectory segments inside 
and outside of the planetary spheres of influence. The trajectory legs 
are constrained to match in position and time- but are initially mis- 
matched in velocity. An iteration scheme is developed to drive this 
mis -match to zero. As the second step, N-body perturbed trajectories 
are calculated which have the same end conditions in position and time 
as the conic legs in the previous step hut have slight offsets in 
initial and final velocities. The same iteration scheme utilized in the 
first step is employed to match these perturbed segments in velocity 
as well as position and time. Finally, the accuracy of each trajectory 
leg is checked by numerical integration. 

Three examples are considered m detail. They are. 

1) a dual planet reconnaissance trajectory 
(Earth-Venus -Mars -Earth 

2) Grand tour trajectory (Earth- Jupiter-Saturn-Uranus - 
Neptune} 

3} periodic trajectory (a repeating Earth-Venus shuttle 
trajectory) 

Free-fall trajectories are determined for the first two of these ex- 
amples. Comparison with numerically integrated trajectory legs has shorn 
these solutions to be accurate to better than 0.4 m/sec m initial and 
final velocity for heliocentric trajectory legs and better than 0.1 
m/sec for planetocentric legs. No free-fall trajectory was found for 
the third example but a powered trajectory (with a total Av of 220.5 
m/sec) is presented. In general, accuracies comparable with the results 
of the preceding two examples are obtained. 
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Chapter 1 
Introduction 

1 . 0 Objectives of the Thesis 

This thesis describes a technique for the determination of accu- 
rate reference trajectories for multiple swingby interplanetary tra- 
jectories. The main objectives of the research are the following: 

1) To develop a basically analytic technique for the deter- 
mination of multiple swingby reference trajectories which 
will converge rapidly from a wide range of initial guesses 
to a solution with a high level of accuracy. 

2) To provide a means of specifying a multiple swingby tra- 
jectory with uniform accuracy along its entire length by 
providing a sequence of guidance aiming points spaced 
along the trajectory 

3) To provide a simple, accurate and economical means for 
performing detailed mission analysis for multiple s\'fingby 
trajectories . 

4) To demonstrate the feasibility, accuracy, and generality 
of the technique by its application to three examples; 

a dual planet reconnaissance trajectory, a Gi'and Tour 
trajectory, and a periodic trajectory segment. 

1 . 1 Definition of the Problem 

The determination of space trajectories is usually posed as a two- 
point boundary value problem. The initial and final position vectors 
and the time of flight between them are given along with the equations 
of motion for the trajectory. The calculation of the reference tra- 
jectory which satisfies these conditions is the targeting problem. This 
thesis deals with the problem of targeting for trajectories charaterized 
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by one or more close planetary encounters between their launch and 
arrival points. 

1 . 2 Existing Targeting Techniques and Their Application to Multiple 

Swingby Trajectories 

The use of 'multiple swingby trajectories to substantially reduce 
the launch energy and flight time for a number of highly interesting 
missions has long been recognized [1,2], Each close planetary encounter 
provides an opportunity to alter the energy of the trajectory with re- 
spect to the sun by use of the planetary gravitational field. In effect, 
the spacecraft exchanges energy with the planet. This ability to make 
major heliocentric velocity changes along the trajectory without fuel 
expenditure allows considerable flexibility in mission planning. Ex- 
amples of some of the missions which have proposed are 

1) Earth -Mars -Earth [3,4] 

2) Earth- Venus -Earth [3] 

3) Deep Space, Solar Probe and Out-of -Ecliptic [5,6] 

4} Earth-Venus -Mercury [7,8,9] 

5) Earth-Venus -Mars -Earth [3,10,11,12,13,14,15] 

6) Outer Planets Missions [16,17,18,19,20,21] 

7) Earth-Venus and Earth-Mars Periodic Orbits [22,23,24,25] 

The majority of these studies have been concerned primarily with 
preliminary mission planning and guidance requirements studies using 
simplified models for targeting. To the author's knowledge, the only 
multiple swingby mission for which precision reference trajectories 
have been generated is the Earth-Venus -Mercury flight [7,8]. 

Present targeting techniques for multiple swingby interplanetary 
trajectories fall into two general classes. The first of these uses an 
approximate model for the mission consisting of a sequence of helio- 
centric conic arcs running from the center of one massless planet to 
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the next. Thus, the trajectory is determined by giving the order in 
which the specific planets are encountered along with the launch, 
arrival, and encounter dates. The effects of the planets on the tra- 
jectory are approximated as impulsive changes m velocity with respect 
to the Sun at each planetary encounter. Using this model, a search is 
made over a range of departure, arrival, and intermediate encounter 
dates to determine combinations which yield trajectories which are 
matched in hyperbolic excess velocity relative to the planet at each 
intermediate encounter and which are physically realizable in the sense 
that they do not require the trajectory to pass beneath the surface 
of any planet. This search may be carried out exhaustively to deter- 
mine all possible swmgby trajectories within the range of dates 
specified, [1,19] or may use an iterative technique to converge on a 
single set of dates [22] . The advantage of this technique lies m the 
speed with which each trajectory may be calculated.. Since the model 
assumes the trajectory to be a sequence of two-body legs, each may be 
determined as the solution to Lambert's Problem. The disadvantages of 
this technique are 

1) The large number of trajectories which must be generated. 
For an exhaustive search procedure, large numbers of 

date combinations must be examined. An iterative technique 
mitigates this difficulty but may not provide all possible 
solutions . 

2) The inaccuracies of the model. Both numerical [7,8,28] 
and analytic [27] studies have indicated that while the 
above model is acceptable for preliminary studies, it 
does not have sufficient accuracy for precise trajectory 
prediction for close planetary encounters. 

In general, approximate targeting schemes presently employed for mul- 
tiple swingby trajectories are most useful foT preliminary mission 
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studies and for the generation of initial conditions for more accurate 

targeting techniques. 

1 

The other class of targeting procedures which have been applied 
to multiple swingby trajectories utilize numerical integration tech- 
niques to generate precision reference trajectories. An example of 
this procedure as applied to an Earth-Venus -Mercury trajectory in [8] 
is as follows, 

1} Initialize the launch conditions at Earth and the aiming 
point at Venus with the conic values from an approximate 
targeting technique. 

2) Search over the injection conditions at Earth until a 
numerically integrated trajectory hits the desired 

aiming point at Venus. 

1 

3) Continue the converged case from 2) on to Mercury and 
note the resulting miss of the desired target point there. 

4) Perturb the aiming point at Venus and repeat steps 2) and 

3). 

5) From the results of step 4), construct partials of the 
miss at Mercury with respect to changes in the aiming 
ppint at Venus. 

6) Compute and apply differential corrections to the aiming 
point at Venus. 

7) Repeat steps 2) 3) , and 6) until convergence at Mercury 
is obtained. 

The average running time for a convergence criterion of + 1000 km at 
Mercury was about 4S minutes on the IBM 7094. A variation on this tech- 
nique employs a.many-body state transition matrix obtained by the 
numerical integration of the variational equations to determine the 
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differential corrections in the targeting process. This latter method 
has been used [26] for the targeting of single leg trajectories (e.g. 
Earth-Venus or Earth-Mars) but to the author's knowledge has not been 
applied successfully to multiple swingby trajectories. 

The advantage of the numerical integration technique is that it 
gives a completely defined accurate reference trajectory. All signif- 
icant disturbing forces may be included to the degree of precision 
available on the computer used. The disadvantages of the technique are 
the large amount of time consumed by the repeated numerical integration 
of the trajectory legs and the question of its feasibility for tra- 
jectories involving more than one intermediate swingby. This latter 
difficulty arises from the strong sensitivity of the trajectory to 
small changes in swingby conditions for planetary encounters earlier 
m the trajectory. Thus, as more swingbys are added to the trajectory, 
the accuracy requirement for the determination of the earlier swingbys 
increases greatly. For the same reason, the linearity region for the 
differential correction process shrinks. Both of these reasons lead 
to a large increase in the number of numerical integrations of tra- 
jectory legs needed. This difficulty did not arise m the approximate 
targeting schemes since the trajectory was modeled as a set of shorter 
arcs to be matched dynamically at a number of intermediate points 
rather than as a single arc determined entirely by its initial con- 
ditions . 

The targeting procedure developed in this thesis attempts to 
combine the advantages of both the approximate and the numerical inte- 
gration techniques while minimizing the disadvantages of both. 

1 ■ 3 Synopsis of Thesis 

In Chapter 2, two patched conic models and their application to 
multiple swingby trajectories are described. The first corresponds to 
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the approximate model described in the preceding section. The second 
(the advanced patched conic model) adds planetocentric conic legs 
between the heliocentric conic legs to describe the swingby maneuver 
more completely and accurately. The conic legs are constrained initial- 
ly to match only in position and time. Then, an iterative process is 
employed to vary the matching points until the legs also match in 
velocity. 

Chapter 3 describes a basically analytic method for computing the 
perturbations of conic legs due to the disturbing accelerations of 
other bodies. A technique is developed for calculating the initial and 
final velocity offsets for each conic leg needed to produce a perturbed 
trajectory having the same initial and final conditions in position 
and time as the unperturbed conic reference leg. 

Chapter 4 deals with the iterative techniques of matching the 
individual trajectory legs (either perturbed or unperturbed) in veloc- 
ity as well as position and time. Both first-order and second-order 
techniques are developed. 

Chapter S presents numerical results for a dual planet recon- 
naissance trajectory. The reference traj'ectory is specified by the 
position, velocity, and time at the sphere of influence entry and 
exit points for the launch, arrival, and swingby planets. Comparison 
with numerically integrated traj'ectory legs indicates that the an- 
alytically calculated trajectory legs match to within a total error 
in velocity of 0.2263 m/sec. 

Chapter 6 presents the same results for a Grand Tour trajectory 
example. Here the trajectory segments were matched analytically to 
within a total error of 2.652 m/sec. 

Chapter 7 discusses a segment of a periodic trajectory that 
shuttles between Earth and Venus. No free-fall trajectory was found 
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for this example but a powered trajectory requiring a total impulse of 
220. S34 m/sec over the 3.6 year segment considered was determined. The 
special nature of this trajectory resulted m less accurate predictions 
by the analytic technique with the total error amounting to 38.950 
m/sec. 

Chapter 8 summarizes the thesis and its contributions. Several 
applications of the techniques developed are suggested for further 
research. 
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Chapter 2 

PATCHED CONIC ANALYSIS 

2 . 0 Chapter Summary 

The application of two patched conic models to multiple swingby 
trajectory analysis is described. The first of these, the simple 
patched conic model, consists of a sequence of heliocentric conic arcs 
matched in relative velocity magnitude at each planetary encounter. 

It is found to be most suitable for preliminary mission analysis. The 
second, the advanced patched conic model, considers the trajectory to 
be approximated by a series of alternating heliocentric and planeto- 
centric conic arcs matched in position, velocity and time at the entry 
and exit points of the sphere of influence of each planet encountered. 
It is found to be a more useful model for reference trajectory calcu- 
lations. The computational details of the advanced patched conic model 
are examined m depth and limitations on its accuracy considered. 

2 . 1 The Simple Patched Conic Model 

The simple patched conic model has been successfully employed for 
a number of preliminary trajectory and mission analysis studies. 
Examples of its use for multiple swingby missions may be found m 
[22, 7, 10, 16, 18, 6] and many others. The model consists of a se- 
quence of heliocentric conic arcs matched in magnitude of velocity 
relative to the planet at each planetary encounter. An illustration 
of one such trajectory is given in Figure 2.1. The steps followed for 
a trajectory with N planetary encounters (launch, N-2 intermediate 
swmgbys, and arrival) are: 

1) Specify the launch date , the arrival date tjg, and the 
N-2 intermediate encounter dates t2> tj,..., t N-i 

2) At each date t^, calculate ■ the position Tp ^ and the 
velocity v p ^ of the planet encountered. 
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3) For each date (except for k=N) calculate a helio- 
centric conic arc running from r p ^ to r p with a 
time of flight T=t^ +1 -t^ (See Appendix B for the method 
of calculating these arcs) . Each arc will have associ- 
ated with it a planetary departure velocity v p ^ and 

a planetary arrival velocity v A both of which are 

measured relative to the sun. 

4) For each intermediate date t^ (k=2 , 3 , . . .N-l) calculate 
the incoming and outgoing velocities relative to the 
planet encountered. 


-1 ,1c — A,k 


' Ip,k 


( 2 . 1 ) 


— 0 ,k ~ — D,k 


- v 




For a free-fall trajectory to be dynamically possible, 
the magnitudes of these velocities relative to the planet 
must be equal at each encounter. Using some convenient 
iteration scheme (see [2Z] for an example) , the inter- 
mediate encounter dates are varied and steps 2 to 4 re- 
peated until this condition is satisfied. 

Note that while the incoming and outgoing velocities (v^ ^ and 

— 0,k^ relative to the planet are equal in magnitude, the arrival and 
departure velocities (v A ^ and v p ^) relative to the sun usually differ 
in both magnitude and direction. The model considers the planetary 
swingby to be equivalent to an instantaneous velocity change of 

A Ik = lD,k - lA,k (2 * 2) 
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relative to the sun applied at the time o£ encounter. 


Once the simple patched conic trajectory has been determined, 
some information on the planetary encounter phases may be obtained. 
Using the incoming and outgoing velocities (v^ ^ and relative 

to the planet as approximations to the asymptotic velocity vectors, 
the constants for a planetocentric hyperbola may be determined and the 
relevant parameters for the swmgby calculated. The accuracy of this 
approximation is studied m [27] . There it is shown that the approx- 
imate swingby parameters differ from their time values by terms of 
order e, the planet-to-sun mass ratio. 

The advantages of the simple patched conic model are its simplic- 
ity, ease of implementation, and speed of computation. A large number 
of trajectory alternatives may be explored with a relatively small 
investment m computer time. Thus this model is well suited for prelim- 
inary mission analysis. The basic limitations of the simple patched 
conic model are: 

1) The heliocentric conic arcs and the planetocentric hyper- 
bolas are matched only approximately. The model does not 
provide a continuous or highly accurate description of 
motion in the vicinity of a planetary encounter. 

2) The effects of the planetary encounter are approximated 
as an impulse rather than considered to act over a region 
in space and time. 

3) All trajectory segments are considered to be conic arcs. 
The effects of all perturbation other than the close 
planetary encounters are ignored completely. 

2 . 2 Advanced Patched Conic Model 

To eliminate some of the inaccuracies and assumptions of the 
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simple patched conic model and to lay the groundwork for a later per- 
turbation analysis, a more advanced patched conic model is necessary . This 
model consists of a sequence of alternating heliocentric and planeto- 
centric conic arcs constrained to match m position and time at the 
sphere of influence (SOI) entry and exit points of each planet en- 
countered. An illustration Is given in Figure 2.2. These entry and 
exit points are chosen initially from the solution for the simple 
patched conic model and usually result in velocity mis -matches between 
the conic arcs. An iterative procedure is necessary to drive this mis- 
match to zero. The SOI used is defined in Appendix C. It is somewhat 
larger than the Laplace SOI commonly used. The application of this 

model to a trajectory with N planetary encounters (launch, N-2 inter- 
mediate swingbys , and arrival) is as follows: 

1) At each intermediate encounter, specify the entry and 
exit points on the planetary SOI. At 'the launch planet, 
specify the exit point and at the arrival point specify 
the entry point. An entry point on the SOI is given by 
its azimuth, elevation, and time of passage. An exit point 
on the SOI is given by the increments in azimuth, eleva- 
tion and time from the corresponding entry point on the 
same SOI. The one exception to this is the exit point 

at the launch planet, which is specified in the same 
way as an entry point. The radius of the SOI is assumed 
to be a constant for each planet. 

2) For each time t^ of entry or exit through the SOI calcu- 
late : 

a) the position Tp ^ and velocity Vp of the planet 
(this is discussed in Appendix G) . 
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b) the cartesian coordinates of the entry or exit 
point with respect to the planet encountered. 


3) For each point k (except the last) calculate a conic arc 
from point k to point k+1 . For even k (entry points), 
this will be planetocentric arc which will always be a 
hyperbola with a central angle greater than 180°. This 
arc will run from the entry point to the exit point of 

a single planet's SOI. For odd k (exit points), the arc 
will be heliocentric and may be either an ellipse or a 
hyperbola. It will run from the SOI exit point of one 
planet to the SOI entrance point of the next. 

4) At each intermediate entry or exit point, calculate the 
difference in velocity between the heliocentric and 
planetocentric arcs. In Chapter 4, an iterative technique 
for varying the position and time of these points in 
order to eliminate the velocity mis-match is described. 

After steps (1) - (4) have been repeated until convergence is. 
achieved, the result is a series of conic arcs continuous in position, 
velocity and time at all points. Discontinuities in acceleration occur 
at the SOI entry and exit points since a planet's gravitational field 
is ignored outside of its SOI while the solar perturbing forces are 
neglected inside of an SOI. Several points to be noted about this 
model are: 

1) For a trajectory with N planetary encounters, there are 
2N-2 matching (entry or exit) points specified in posi- 
tion and time along the trajectory. Odd-numbered matching 
points are SOI exit points while even-numbered ones are 
SOI entry points. 
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2) The different means used to specify entry and exit points 
provides some separation in the effects of varying these 
points. The two heliocentric arcs touching the planet's 
SOI are affected primarily by changes in the entry point. 
The heliocentric arc leaving the planet is, also affected 
by exit point changes but these effects are usually much 
smaller than those due to the entry point changes. The 
planetocentric arcs are affected primarily by exit point 
changes . 

3) This model considers the effects of an planetary encoun- 
ter to be distributed over a region in both time and 
space. This is a more accurate approximation to the ac- 
tual interaction than is provided by the simple patched 
conic model. 

4) A continuous description of the motioh along every phase 
of the trajectory is given. This also provides a basis 
for the perturbation analysis to be described in Chapter 

3. 

A detailed discussion of the computations involved in several of 
the steps in the advanced patched conic model is given in the next 
section. 

2 . 3 Computational Details 
2 . 31 Entry and Exit Point Coordinates 

The state vector for an entry or exit point is given by 


*k 



k=even 
or k=l 


(entry point) (2.3) 

(launch point) 
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k=odd and ^ 1 


(exit point] 


(2.4] 


£k = 


&0 k 

A<J>j 

At, 


where 9^ = azimuth of point k 

= elevation of point k 
t k = time of passage through point k 
A6j, = difference xn azimuth between points k and k-1 
A^ = difference in elevation between points k and k-1 
At^ = difference in passage time between points k and k-1 


These coordinates are illustrated in Figure 2.3. 

The planetocentric cartesian coordinates of an entry or exit point are 
calculated using 


cos a cos P 


I 

a 

L* 

sin a cos P 

= 

,k 

sin 0 


_ z D,k_ 


where r . = radius of SOI for planet i 

S 9 1 

a = 0^ k=even or k=l 

= e ]c _ 1 + A6^ k=odd and k^l 
g = (j>^ k=even or k=l 

= ^ 1 + Atfi-^ k=odd and k? 5 ! 


(2.5] 


2 , 32 Calculation of Conic Arcs 

In [3], it is shown that given their initial (r^, t^] , and final 
(r 2 , t 2 ) positions and times, it is possible using Lambert's theorem 
to calculate a two-body conic trajectory with initial velocity v^ and 
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final velocity v 2 connecting any two points. This procedure is 
described in detail in Appendix B. 

2 . 321 Planetocentric Arcs 

These arcs run from the entry point on a planet’s SOI to the 
corresponding exit point. They are always hyperbolic with respect to 
the planet and traverse a central angle greater than 180° but never 
make a complete revolution. The initial and final points are 



(k=even) (2.6) 

—2 = — Djk^l H = *k + At k*l 

The initial and final velocities are stored as 


-H,k = -1 


(2.7) 


— H,k+1 = —2 


2 .322 Heliocentric Arcs 

The heliocentric arcs run from the SOI exit point on one planet's 


SOI to the entry point of the next. 

The initial 

and 

final points are 

*1 = Ip,k + — D,k 

*1 = t k-l + 

At k 

(=t^ for k=l) 

-2 = — P,k+1 + -D,k+1 

t 2 = t k+l 


(k=odd) (2.8) 

In addition it is necessary to 

specify 




1) The number of complete revolutions the arc traverses 
about the central body. This must be given a priori. 
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2) Whether that portion of the arc remaining after the 
complete revolutions have been finished traverses a 
central angle greater or less than 180°. Assuming that 
all heliocentric trajectories have inclinations less 
than 90°, this may be determined using 

Gj_ = sgn [i z *(r 1 xr 2 )] (2.9) 

If (^>0, the central angle is less than 180° while if 
G-^<0 it is greater than 18D° . 

3) Whether the arc is a hyperbola or an ellipse. This is 
determined by comparing the time of flight for the arc 

T * t 2 - t 2 (2.10) 

with the parabolic time of flight (see [3]) 

!p= \ Vi*[ s3/2 ' g i^- c ) 3/2 ] ( 2 • X1 ) 

£ = I 2 ’-l c = 1-1 

s = 2 ( r 1 +r 2 + c) 

p = gravitational parameter for the central body 
(m this case the sun) 

G^= as defined above 

between the same two points. Then, if T>Tp, the arc is 
an ellipse while if T<T p , it is a hyperbola. If T=Tp, the 
arc is a parabola. 
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The initial and final velocities for the heliocentric arcs are 


stored as 


i E> k = h 

— E ,k+l = —2 


( 2 . 12 ) 


2 ,53 Calculation of Cost Function 

At each entry or exit point along the trajectory (excluding the 
initial and final points) the velocity along both a heliocentric arc 
and a planetocentric arc have been calculated. In general, these 
velocities will not be consistent but instead will be mis-matched by 
an amount 


A ^k = — E,k 


^H,k ‘ -P,k 


(2.13) 


where — E k = velocit y relative to the sun along the heliocentric 

arc at point k 

“H , k = velocity relative to the planet along the planeto- 
centric arc at point k 

Vp k = velocity of the planet relative to the sun at time 


A scalar cost function J for the total velocity mis-match along 
the trajectory may then be calculated using 


2N-3 _ 

J - l AV Av^ (2.14) 

k— 2 
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This expression is positive definite and goes to zero only when 
velocities are matched along the entire trajectory. 

2 . 4 Accuracy of the Model 

The basic limitation on the accuracy of this model lies m the 
fact that it assumes each trajectory segment to be a pure two-body arc. 
The effects of direct planetary gravitational attractions are ignored 
outside of the planet's SOI while the effects of solar perturbing 
forces (due to the gradient of the sun's gravitational field) are neg- 
lected inside of a planet's SOI. Within these assumptions, the calcu- 
lation of the trajectory segments is an exact solution to the non- 
linear two-body orbit determination problem. The accuracy of the so- 
lution is limited only by the computational round-off errors of the 
method used. 
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Chapter 3 

Perturbed Conic Analysis 

5 . 0 Chapter Summary 

Two approaches to computing the perturbed conic trajectory seg- 
ments are considered. The first computes the perturbations for a tra- 
jectory running from the initial time to the final time. Perturbations 
due to disturbing accelerations near the initial time are found to 
grow to unacceptable levels near the final time. The second method 
starts at the trajectory mid-point and computes perturbations from 
there to both the initial and final times . The accuracy for this meth- 
od is adequate. The details of evaluating the perturbations by quad- 
rature are described. The method for determining the velocity offsets 
for a perturbed trajectory passing through the same initial and final 
position and time as the two-body trajectory is given. The source of 
error in the calculations are discussed. 

3 . 1 Description of Approach 

The object of this chapter is to take into account the fact that 
the segments of the true multi-swingby trajectory are only approxi- 
mately two-body orbits. To do this 3 a set of perturbed conic trajectory 
segments (corresponding to the segments of the advanced patched conic 
model) are calculated. The heliocentric arcs take into account the 
disturbing effects of the attraction of the planets while the planeto- 
centric arcs are affected by the disturbing forces due to the sun. The 
perturbed conic trajectory segments are constrained to match the same 
initial and final conditions in position and time as the corresponding 
segments m the advanced patched conic model but are offset m initial 
and final velocity. It is the calculation of these velocity offsets 
that is the mam concern of this chapter. 

Two approaches to the calculation of the velocity offsets were 
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tried. The first of these is similar to the implicit velocity offset 
technique used in space guidance [28]. This procedure is illustrated 
in Figure 3.1. Its steps for a single trajectory segment are described 
below. 


1) Assume that the state 


— 0 


loctr 

.Soft). 


(3.1) 


along the advanced patched conic trajectory segment found 
m the preceding chapter may be expressed as a function 
of time. Starting with the same initial conditions as 
the two -body trajectory 


xCt-^ = 


(3.2) 


calculate the perturbed trajectory x(t) using linear 
perturbation theory. 

2) At the final time, calculate the position and velocity 
differences between the two-body and perturbed trajec- 
tories . 


6x(t 2 ) 


6l(t 2 ) 

SXCV. 


= x(t 2 ) - x Q (t 2 ) 


(3.3) 


3) Using linear perturbation theory, calculate the velocity 
offset needed at the initial time to reduce the final 
time position offset value to zero. 
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Before After 



0 



0 

6x(t,) = 



Sx(t-) = 



_0 



6v(ti) 

Sx(t 2 ) = 

6r(t 2 ) 

_6vCt 2 )_ 

6x(t 2 ) = 

0 

dv(t 2 ) 


This process also leads to a new velocity offset value 5v 2 (t 2 ) 
at the final time. The above procedure proved to be highly inaccurate. 
All the trajectory segments have in common the characteristic of moving 
from a region of strong perturbing forces [near the sphere of influence 
boundary) into a region of weak perturbing forces (far from any planet 
for the heliocentric legs and deep inside the sphere of influence for 
the planetocentric legs) and then back into a region with strong per- 
turbing forces. The effects of the initial strong perturbations grow 
rapidly along the trajectory and lead to very large deviations in 
position and velocity at the endpoint. These deviations are usually 
outside the linear range of the perturbation theory used in steps 1 
and 3 above , thus leading to unacceptable errors in the calculation 
of the velocity offsets. 

The reason for the failure of this first approach led to the use 
of the much more successful second approach. Referring to Figure 3.2, the 
steps for this method for a single trajectory segment are given below. 

1) Assume that the state Xg(t) along the advanced patched 
conic trajectory segment is a known function of time. 
Starting at the mid-point 

tj^ = 2 (^2 + (3.5) 

of the two-body segment with the same state as the 
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conic trajectory, calculate perturbed trajectories 
forward to the final time and backward to the initial 
time (See figure 3.2a] from that point. 

2) At the initial and final points, calculate the offsets 
in position and velocity between the two-body trajectory 
and the perturbed trajectory. 


dxCt^ 


Sx(t 2 ) 


drftj.) 

<$v(t,) 
- X - 

'<Sr(t 2 )‘ 

<Sv(t 2 ) 


x(t x ) - XgCt-L) 


x(t 2 ) - 


(3.6) 


(3.7) 


3) Using linear perturbation theory, calculate the offsets 
in both position and velocity at the mid-point time t^j 
needed to reduce the initial and final position offsets 
to zero . 


Before 


After 


6x (t^) = 


<Sx(t M ) = 


6x(t 2 ) = 


6r(t 1 ) 

6 ^V| 


<sr(t 2 ) 

dvCt 2 )_ 


6x(t 1 ) = 


= 


6x(t 2 ) = 


0 

6v( tl ) 

«i(V 

L 5v(t M )j 

0 

6v(t 2 ) 


(3.8) 


This new trajectory is shown in Figure 3.2b. 
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This second approach provided the accuracy needed for the calcu- 
lation of the initial and final velocity offsets (<5v(t^) and <Sv(t2)). 
Since the perturbed trajectories calculated in step 1 always run from 
regions of weak perturbations (the mid-section of the trajectory) to 
regions of strong perturbations (the end-points), large deviations 
(due to the accumulated effects of the strong perturbing forces) do 
not have the time to grow. Also, since the effects of the strong per- 
turbing forces depend on the time spent in the vicinity of the sphere 
of influence boundaries at the end-points rather than on the time of 
flight of the trajectory, the size of the position and velocity offsets 
(6x(t^) and 6x(t 2 )) are not influenced heavily by the length of the 
trajectory. The computational details of the second approach are de- 
scribed in the next section. 

Once the velocity offsets have been determined, a new cost func- 
tion taking them into account is contructed. Then, the iterative 
procedure employed fox the advanced patched conic model is used to 
match the perturbed conic segments in velocity as well as position and 
time at the entry and exit points. 

3 . 2 Computational Details 

5 . 21 Calculation of Perturbed Trajectory 

It may be shown (3,29,30,31] that the solution (using linear per- 
turbation theory) for the deviations between the perturbed and two-body 
trajectories is given by 

t 

<5x(t) = $ 0 (t,t..) 6x(t i ) + f <f 0 (t,r) ^(t) dr (3.9) 

*i 

where 
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rsr(t)1 




L<Sv(t)J 


= deviation from two-body orbit 




- 


0 

5d^) 


= state transition matrix for the two- 
body orbit between t^ and t. (See Ap 
pendix D for a description of the 

properties of this matrix.) 

= disturbing vector 


a^(t) - = disturbing acceleration evaluated as a 

function of time along the two -body 
traj ectory. 


All the quantities on the right-hand side of (3.9) are known functions 
of time evaluated along the two -body reference trajectory. Since the 
perturbed traj'ectory calculations start with 


x(t M ) ~ x Q (t M ) 


it can be seen that 


«x(t M ) «■ 


0 

0 


and that (3.9) becomes 


fixCti) 



V t i, T )Io COdT 


(3.10) 


(3.11) 


(3.12) 
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for the integration to the initial point and 


i$xCt 2 ) 



Ct 2 ,T)f 0 Ct)dT 


(3.13) 


for the integration to the final point. Since the integrands in (3.12) 
and (3.13) are known functions of x, the integrals may be evaluated by 
quadrature (See Appendix D for a discussion of the technique used) . 

3 .22 Calculation of Disturbing Accelerations 

The disturbing acceleration due to body P^ on the motion of P 2 
with respect to P.^ (see Figure 3.3) is given by 


where 



JU = position vector from P^ to P^ 


(3.14) 




position vector from P^ to P 2 


Uj = gravitational parameter for P^ 


Numerical difficulties may arise in the use of (3.14) since A and d^ 
are often nearly equal and opposite vectors. These difficulties may be 
alleviated using a technique developed in [3] . Write 



(3. IS) 


where 
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Now, write 


r = position vector from P^ to T?^ 


h * 


(S' 1 )' wtq - 


where 


*j ‘ ij (i, ' 2 “ s “j) 


w Cq j ) = Ci + q.j ] 3/2 - i 


a = angle between r and &. 
3 6 - -3 


To evaluate W (q^ ) , re-write (3.18) as 


(l + Q-i J _ 1 

W(q ) = — r H2 

3 (l+q^^+l 


or 


3 + 3q +qf 

W(q.= ) = q. zh — 

J J (l+q 3 0 ' + 1 


Thus, substituting (3 . 19) and (3 . 16) into (3.15) yields 




.j ' • ^ [i * wf V k 


(3.16) 

(3.17) 

(3.18) 


(3.19) 

(3.2Q) 
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The total disturbing acceleration is the sum of the individual 
contributions . 




N 

= l fid 
3=1 


(3.21) 


For the heliocentric legs, all significant planetary disturbing accel- 
erations are included in the calculation of (3.21). For the planeto- 
centric legs, only the disturbing acceleration due to the sun is con- 
sidered. 


5 .23 Calculation of Velocity Offsets 


As shown in Appendix D, the state transition matrix may be par- 
titioned into four sub -matrices . 


where 


= 


AoCt.tf) 


A Q (t,t i ) 


c 0 (t,t.) 


3r(t) 


B 0 (t,t.) 




(3.22) 


B 0 Ct,ti) - 


3r(t) 

LSvUTJ 


Cjj(t,t i ) = 


3v(t) 
3r(t i )_ 


D 0 (t 


favct) 

*V ~ b^o. 


Having computed the initial and final perturbations (6x(t 1 ) and 6x(t 2 )) 
corresponding to 6x(t M ) = f° rm the correction matrix 
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H 0 ^ t l» t 2 ,t M^ 


A 0 (t l»V 

A 0 ( ' t 2 ,t M- ) 




(3.23) 


The new value for the deviation of the state at the mid-point <$x(t M ) 
is calculated to reduce the position deviations at the initial and 
final points to zero. From (3.9), (3.12), (3.13), (3.22), and (3.23), 
it can be seen that 


Thus , 


’firCV 


£ 

6r(t 2 )_ 


£ 




6r(ti) 

6r(t 2 ) 


(3.24) 


h q Si- 1 


6r (ti) 
6r (t 2 ) 


(3.25) 


The new values for the state deviations at any point may be cal- 
culated using (3.9). 


t 

6x(t) = $ 0 (t,t M ) <5x(t M ) + J 4> c (t,T)f Q (t)dt (3.26) 

t M 


Specifically, the offsets at the initial and final points are given by 


dx(ti) 


£ 

fivCtj) 




(3.27) 
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6x(t: 2 ) = 


0 

Sv(t 2 ) 

3.24 Calculation of Cost Function 


= vw 4cv + <5x(t 2 


For a heliocentric arc running from t^ = t^ to t 2 
the velocity offsets as 


ov 


-E,k = *£<* 1> 


(k = odd) 


Sv 


i El k + i = 


For a planetocentric arc running from t^ = t^ to t 2 
the velocity offsets as 


61^^ = 6v( tl ) 

%,k + l = 


(k = even) 


Then, let 


A ^k = Xs.k + E ,k ' ^P,k - ^H,k 


6 XH,k 


and 


2N-3 


J = J„ A Xk T A ^k 


k=2 


) (3.28) 


Vi* store 


(3.29) 


= Vi* store 


(3.30) 


(3.31) 


(3.32) 
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3.3 Sources of Error 


The errors associated with the calculation of the perturbed tra- 
jectories have as their source 

1) Computational errors, m such areas as matrix inversion, 
calculation of the two body orbits and state transition 
matrices, evaluation of the perturbation integrals by 
quadrature, etc. These errors may be reduced to any level 
desired by increasing the precision of the calculations, 
reducing step-size fox quadrature methods, and by in- 
creasing the accuracy level required for the termination 
of iterative solutions to transcendental equations. The 
ultimate level of accuracy due to errors of the above 
nature is limited only by the precision available on the 
computer used. 

2) Errors associated with the approximations involved m the 
use of the trajectory model. Referring to (3. 9], it can 
be seen that these errors occur because of the 

i) evaluation of the disturbing acceleration on the 
two -body reference trajectory rather than on the 
perturbed trajectory. 

ii) use of the two-body state transition matrix 
calculated for the reference trajectory rather than 
the many-body state transition matrix §(t,t^) eval- 
uated along the perturbed trajectory. 

lii) use of linear perturbation theory. 
iv) neglect of smaller disturbing forces. 

The error sources listed under 2) above are the dominant factors. 
Examining each of these sources in detail their importance can be es- 
timated. 
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i) From (3.14), the disturbing acceleration is 



The gradient of this acceleration is given by 



where is considered a constant. 


(3.14) 


(3.33) 


The change in a, - due to a small shift from the reference trajectory 
~d , 3 

is given by 


6 Sd , 3 = G ^j) % 

= G(d ) 6r (3.34) 

since r = & + implies 6r = . From (3.33) and (3.34), it can be 

seen that the magnitude of 5a. is roughly 

> J 



The largest deviations from the reference trajectory for either planeto- 
centric or heliocentric legs occurs at the sphere of influence (SOI) 
boundary. From the definition of the SOI (see Appendix C) , the dis- 
turbing acceleration due to the sun equals the primary acceleration due 
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to the planet on this surface. On the heliocentric legs, the disturbing 
acceleration at the SOI boundary is largely due to the primary accel- 
eration of the planet while on the planetocentric legs, the disturbing 
acceleration at the SOI boundary is due largely to the perturbing ac- 
celeration of the sun. From (C.9) , the disturbing acceleration at the 
SOI boundary may be written approximately as 

' (?) r ' (if r ) (3 ' 36) 

where 


y = gravitational parameter of the planet 
Pj = gravitational parameter of the sun 
Since, at the SOI, r<<£, (3.36) may be written as 


a 


d.j 



(3.37) 


Then, from (3.35) 


(3.38) 

a d,0 T 

at the SOI boundary. Typical maximum values for this ratio are 0.06 
for inner planets and 0.01 for outer planets. 

ii) Numerical studies comparing an analytic two-body state transi- 
tion matrix with a many-body state transition determined by numerical 
integration showed that terms in both matrices remained equal to 
within a few percent for both heliocentric and planetocentric legs 
when the enlarged sphere of influence (see Appendix C) was used. When 
the Laplace SOI was used, large differences (often over 100%) occurred 
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in terms m the C and D matrices for the heliocentric legs. 

ni) Numerical studies indicate that the perturbations ecountered 
in all the steps are small enough for the linear theory to remain valid. 
This is also shown by Slater and Stern in [28] . 

xv) The perturbing forces neglected m this model are those due 
to other planets during planetocentric legs, those due to oblateness 
and other higher-order terms in the gravitational field of the sun and 
planets and those due to non-gravitational effects such as drag, radi- 
ation pressure, etc. They are considered small compared to the forces 
included. 

Numerical values for the accuracy of the perturbed trajectory 
calculations may be found m Chapters 5-7 and Appendix D. In general, 
the perturbed trajectory appears to eliminate about 98%-99? of the 
error between the conic and N-body trajectory segments. The largest 
error source is usually the difference between the two-body and N-body 
state transition matrices. This result is also indicated m [28]. 
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Figure 3.1 Offset Calculation f^First Method 



Figure 3.2 Offset Calculation [second method! 
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Figure 3.3 Disturbing Acceleration Geometry 
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Chapter 4 


Trajectory Segment Matching Procedure 

4 . 0 Chapter Summary 

The problem of matching the trajectory segments in velocity as 
well as position and time is formulated as a parameter optimization 
exercise. The details of the calculation of the gradient of the cost 
function for the velocity mis-match are described. A number of first- 
order techniques (steepest decent, modified steepest descent, con- 
jugate gradient, and acceleration steps) are applied to the problem. 

A second-order technique (generalized Newton-Rapheson) is also dis- 
cussed and applied. The behavioT of the different techniques is des- 
cribed and the best are selected. The application of inequality con- 
straints on the distance of closest approach to each planet is de- 
tailed, The application of the trajectory segment matching procedure 
for the multiple swingby analysis is outlined. 

4 . 1 Description of the Problem 

The problem of minimizing the velocity mis-match at the trajectory 
entry and exit points may be formulated as a parameter optimization 
problem. For a trajectory with N planetary encounters (launch, N-2 
intermediate swingbys, and arrival) , define an expanded state vector 
and its variation as 


—2 



—3 


6x 3 

• 


* 

• 

6s = 

• 

• 


* 

-2N-3 


5 -2N-3 

__ . 


— — 


where Xj, is defined as m (2.3) and (2.4). The launch point x^ and 
arrival point ^2N-2 are cons id-e r ed fixed, so that 
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6 *1 = 2N-2 = 5. 


(4.2) 


Similarly, for the 2N-4 trajectory matching points (the intermediate 
entry and exit points), define a velocity mis -match vector as 


Av„ 


^—3 


2N-3. 


(4.3) 


and a cost function 


2N-3 T 
J = l Av k Av. 
k=2 K K 


(4.4) 


where Av^ may be defined either as m (2.13) or in (3.31) 

The object of the problem is to find the value of which min- 
imizes J and if possible reduces J and u to zero. Both <$£ and u have 
6(N-2) independent components, so that sufficient degrees of freedom 
exist for a solution to be possible. 

To avoid unrealizable trajectories, it is necessary to apply the 
inequality constraint 


r CA * k M 


(4.5) 


where 


r^ = radius of closest approach of trajectory to the planet 


54 



r g = equatorial radius of the planet 

= constant multiplier (nominal value = 1.1) 


at each swmgby. Changes in the value of <$s which violate this inequal- 
ity constraint are not allowed. The details of the constraint appli- 
cation are discussed m Section 4.5 of this chapter. 

4 . 2 Calculation of Cost Function Gradient 

4 . 21 Calculation of Lambert Problem Partials 

During the calculation of the conic trajectory segments for the 
advanced patched conic model, it is also possible to calculate analyt- 
ically the partial derivative matrices 



3v. 


Vl 

' * 

(4.6) 


3V n 


^2—1 

_ _L 

ar 2 

(4.7) 


3v 0 


^1—2 

— Z 

(4.8) 


3v~ 


V 2— 2 

— Z 

= ^ 

(4.9) 


and the partial derivative vectors 


3v 2 3v x 3v 2 
clt^ ’ 3t x ’ 3t 2 ’ 3t 2 


(4.10) 
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where — 1’ t l = initial position, velocity, and time for the 

conic arc 

r 2 , v 2 , t 2 = final position, velocity, and time for the 
conic arc 

The relations necessary for the calculation of these partial derivatives 
are given in Appendix B. They aTe in cartesian coordinates and must be 
converted to spherical coordinates for use m the cost function gra- 
dient. Recall from (2.6) that for a planetocentric conic segment 


£l = — D,k H = t k 

-2 = -D,k+1 t 2 = t k + At k+1 

and from (2.8) that for a heliocentric conic segment 

-1 = k * -D, k = tj- + At k 

(= t fe for k=l) 

—2 ~ — P, k+1 + -D, k+1 t 2 = t k+l 

where, from (2.5), it may be written 


(4.11) 


(4.12) 
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Since motion of the entry and exit points on the sphere of influence 
(SOI) does not affect the planetary position r p ^ 


% 


D ,k 




5 -D,k+l 


(4.14) 


for both planetocentric and heliocentric arcs. The relation between 
variations m cartesian and spherical coordinates is given by 


6r = 
-J 


67 1 

L 6z i. 


- r s,i co *+j 

sin0 ^ 

_r *,i 

sin<J>j cos6^ 

cos4>j 

COS 0 

1 

r s,i cos *j 

COS0j 

* r s,i 

sin^j sin0j 

COS^ 

sme^ 

0 


r 3,l 

cost{> 


simj). 
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a *j 


,6r 


(4. IS) 


Since the entry and exit points may vary only on the SOI, it is 
necessary that 


5r s,i E 0 


(4.16) 


so that (4.15) may be written as 
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-r cos* sin0 - 

s > 1 J J 

r cos* cos6. 

S 7 1 J J 

0 


= 


50j 

5<i> n 


L o J 


„ ”1 


— 

r Sin6 . cos6 0 

^9 ± J J 


de. 

r - sin<j>. sine. 0 

s > 1 j j 


S<t , 

r , cosd). 0 


0 

S,1 


_ 


(4.17) 


Referring to (4.11) and (4.12), it can be seen that the partial 
derivative vectors (4.10) are affected by planetary motion for the 
heliocentric legs but not for planetocentric legs. This effect is 
given by 


d -i !zi + i=i 

' 3t i 3 ^p,k st i 


3Vi 


at^ + tv i-i) ~p,k 


(4.18) 


since 


5l 

9 *P,k 


^—1 

3r7 


7 lll 


when — D k = fixed. 

Similarly, 


dv^ ‘Sv^ 

dt^ jq + Xp, k+ i 


(4.19) 
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dv_2 9 v 2 

3t^ = 9t^ + ^1—2^ — P,k (4.20) 

dv ? 3v ? 

“ Fq + ^2—2^ ^P,k+1 (4,21) 

Note that the relations (4 . 18) - (4 . 21) are used only for the helio- 
centric legs. For the planetocentric legs 


dv. 

^7 



Cx, 3 = 1, 2) 


(4.22) 


If a new state vector y_ (differing from the state vector x) is 
defined as 





(4.23) 


for all j , it is then possible to define the partial derivative 
matrices 


^1 

3 Z.l 


" Qn " R i + 


i^-i J i 


dvi 


(4.24) 


S t2 


Q21 ~ (^?Kl) R 7 + 


2~-l J 2 




(4.25) 
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3v n 


flv “ Qx2 = R 1 + 


3 *1 


1—2 


0 

0 

0 


d — 2 

at^ 


(4.26) 


a I 2 

3^2 


~ ^22 “ ^ 2 - 2 -* R 2 + 


d — 2 


(4.27) 


For each conic arc (either heliocentric or planetocentric) from 
point k to point k+1, the partial derivatives (4 . 24) - (4 .27) are stored 
as 


A k = 

Q 12 

(4.28) 

B k + 1 

= ^22 

(4.29) 

c k - 

Qn 

(4.30) 

D k+l 

= ^21 

(4-31) 


4.22 Calculation of Gradient 

From (2.14), the cost function was given by 


2N-3 _ 

J = I Av, Av v 
k=2 ~ K ~~ K 


(4.32) 


Its gradient with respect to x^ is given by 


^ 2N-3 T 3Av, 

_ 3J _ „ v 1 — k 

z k=2 ~ k 3 -Jl 


(4.33) 
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where 


H = 2,3,4 , .... 2 1-3 
and, from (2.13) 




— E,k ' — P,k 


^H,k 


(4.34) 


From the definition of the state vector x^ given m (2.3) and (2.4), 
the effects on the state vector defined in (4.23) and hence on the 
terms in (4.34) due to a change in Xj, occur at 


Point affected 


£k 


-E ,k 


-P,k 


^H,k 


k-1 


k+1 

k+2 

for an entry point and at 


Point affected 


Xk — E ,k ^P,k 2fl,k 


k-1 


k 

k+1 


for an exit point. The effect on Vp ^ is due solely to the change in 
the time t^ associated with point k. This is given by 
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p 


k 


^P,k 

3 ^k 


'o 

0 

1 - 


■ 0 

0 

1 - 

0 

0 

3 ^p,k 

3t k 

= 

0 

0 

— P > k 

.0 

0 

1 - 


. 0 

0 

1 _ 


where 

-P,k ' 3 -P,k 

r P,k 

U s = gravitational parameter of the sun 
^ « as defined in (4.23) 


(4.35) 


(4.36) 


Thus, using (4.33) the gradient at an entry point may be written as 




( T 
2 < Av 


— k-1 


^E.k-l 


3 *k 


+ A V 


3v -E,k _ a £p,k 


|_ 3 -k 


9x, 


3 Sk 


+ A ^k*l 


3 ^E,k+l _ 3 -P,k+l _ 3 ^h,k+l 


3 *k 


3 *k 


3 — k 


- a ^; 2 


3 —E,k+2 


L 3 — k 


]} 


(4.37) 


Using the relations (4. 28) - (4. 31) , this becomes 


82 



s-k - 2 | AtJj [D k ] . Av k T [b* - p k - c k - D kn ] 



T 

+ A ^k+1 

[- A k - Vi - p k.i * S.i] 

(4.38) 


, . T I 
+ A — k+2 

[vj} 


9 “)> 

9 “* 

. ^ 


9 *k 

9 *k 

9i k+l 



-Z ~ — E ,Z ’ —P,Z 3 ° r , Z 


for an entry point. 

For an exit point, the gradient is given by 


u _ 

^■k 



9 — H , k~ 
3x k 


+ A ^Jl 



(4.39) 


Using (4 . 28) - (4 . 31) , this becomes 


«k 



[-»k] * a ^ t [W’k] 


- **Ji MS' 


(4.40) 


63 



since 


for an exit point 





The gradient vectors are calculated for all entry and exit points 
except the launch and arrival points (these are considered fixed bound- 
ary conditions) and then formed into an enlarged gradient vector 


% 


- s 2N-3 


(4.41) 


4 . 3 First-Order Techniques 

In finding a solution z=ct of an equation f(z)=a, an iterative 
technique which functions such that 


Z k + 1 - “ = q k ( V a ) 


|q k ! < 1 


(4.42) 


is known as a first-order technique. Several first-order techniques 
were employed to minimize the ’cost function J. 

4 . 31 Steepest Descent 

Referring to the notation of (4.1), the method of steepest 
descent prescribes a change 


— k+1 = \ + 6 *k C4 ‘ 43) 

«£ k = - \ £(s k ) (4.44) 
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The multiplier h k is used to set the step size. One method of 
determining it is on the basis of desired change in the cost J. Since 


6J = £ (£ k ) 

= ' h k & (! k ) 


J 

2. 

the value of h, is given by 


(4.45) 


6J 


k . £ T CfLk^k> 


(4.46) 


The value of 6J is chosen to be some fraction of the cost J at the 


state s^. Thus 


Thus 


SJ ~ - Y k J (^. k ) 

V c ^k> 


h k = ~1 


(4.47) 


(4.48) 


£ Cs k )g(s k ) 

The initial value of y k is y^ S 1 • 

The procedure for a single steepest descent step is as follows 

1) At the present state s k , evaluate the gradient g.(£ k ) and 
the multiplier h k< 

2) Take the step given by (4.44). At the new state s k ^, 
evaluate the cost >Hs k+1 ) • 

3) If the new cost JCs^+p) is less than the old cost J (£ k ) , 
accept the step and set 


^k+l ~ ^k 


(4.49) 


4) If the new cost is greater than the old cost, reject the 
step. Set 


Y k +1 “ °’ 5 Y k 


(4.50) 
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and repeat the step . 

This form of the steepest descent technique suffers from two 
shortcomings. The first of these is the method of step-size control. 
While it does insure that each accepted step will reduce the cost 
function, it does not usually take the best possible step. The second, 
and more serious, concerns the nature of the cost function, which is 
much steeper m directions corresponding to changes in the angle com- 
ponents of s_ than in directions corresponding to changes in the time 
components. This is the common "ravine" problem encountered in many 
parameter optimization situations. As depicted in Figure 4.1, it 
results in a zig-zag path yielding little cost reduction for each step. 
The means of alleviating both these shortcomings are dealt with in the 
next section. 

4 .32 Modified Steepest Descent 

4.321 Optimum Step-Size Selection 

To insure an optimum step in the direction specified by (4.44), a 
parabola is fitted to the cost function in that direction and the step 
taken to its minimum. For a parabola given by 

(y-c Q ) = c iO-c 2 ) 2 (4.51) 

the constants Cg, c-^, c 2 may be determined from the values of the 
ordinate y at x=0 and x=l and the slope at x=0. 

*0 • ?<“> V * £ I „ w- 52 ’ 

x = 0 

7 1 = yd) 
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Substituting (4.52) into (4-51) yields three equations 


y 0" C 0 ~ c l c 2 2 

(4.53) 

y l“ c 0 = c 1 d-c 2 ) 2 

(4.54) 

V = - 2c 1 c 2 

(4.55) 


Solving (4 . 53) - (4 . 55) for Cq, c ^, c 2 yields 


c i = y r>W 



(4.56) 

(4.57) 


y 0" C l C 2 


(4.58) 


This technique searches for the minimum of the post function 
along the line £ = £ k +x6s. . The units are such that x=l corresponds 
to the step taken in the preceding section. 

The procedure for taking an optimum step is as follows: 

1) At the present state s^, evaluate the cost ^{s^), the 
gradient £(£ k ) , and the multiplier h^- 

2) Take the step 6£ k given by (4,44), At the new state 
£ k evaluate the cost J(£ k ). 

3) For the parabolic fit, set 

7 0 = J (s k ) (x=0) (4.59) 
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Cx=l) 


(4.60] 


7l = J(£ k ) 

V = 5J " ' Y k J(s k ) (x=0) (4.613 

Solve for the values of Cg, c-^ , c 2 using (4 . 563 - (4 . 58} . 
43 Take the optimum step 

-k+1 a ik + c 2 6 *k ( 4 ' 6Z ) 

and set 

Y k+ 1 = c 2 y k (4.633 

This technique requires one additional evaluation of the cost function 
(in step 23 

4.322 Gradient Weighting Matrix Selection 

From Figure 4.1, it can be seen that the gradient components 
directed across the ’'ravine" will experience a sign change at each 
step while those components along the axis of the "ravine" will remain 
unchanged in sign. Thus, the ravine problem may be alleviated by 
adding a weighting matrix W^. which turns the step direction along 
those gradient components which do not change sign. This modification 
of (4.443 results in the new step 

<*l k = ' h k W k &CiLk) (4.643 

The weighting matrix is a diagonal matrix of weighting coefficients 
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w, 


l,k 


w 


2 ,k 


w- 


3,k 


(4.65) 


W 6N-12,k 


each o£ which corresponds to one component o£ the gradient vector. 
The initial value of the weighting matrix is an identity matrix. 


W x = I (4.6$) 

After each step, each component of the gradient is tested to see if 
it has changed sign. The weighting coefficients corresponding to 
those components which have not changed sign are increased by an 
amount w„. 


w. . = w. . + w 

i,k+l i,k g 


(4.67) 


while the coefficients corresponding to components which have changed 
sign are decreased by a like amount. 


w. , = w. , - w 

i,k+l i,k g 


( 4 . 68 ) 


All weighting coefficients axe constrained to lie in the range 


1 ± w i,k 5 W M 


(4.69) 
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Different values of gain w^ and maximum value w^ are assigned to 
coefficients corresponding to the S, <j>, and t components. 

The net result of this procedure is to bias the direction of the 
step in favor of those components of the gradient whose sign remains 
unchanged. This is exactly what is needed to proceed down the length 
of the "ravine". The technique which incorporates both the optimum 
step-size control and the weighted gradient step will be termed the 
optimum gradient step (OGS) . 

4 .33 Conjugate Gradient Method 

The method of conjugate gradients uses the information gained 
from previous steps taken to gradually construct a set of mutually 
conjugate directions. If the cost function J were an N-dimensional 
quadratic form, a sequence of N one -dimensional minimizations along 
these conjugate directions would locate the minimum of the cost 
function. For general functions, the process is iterative rather 
than convergent m a finite number of steps. The method is detailed 
in [32]. The computer subroutine used is DFMCG taken from [33]. 

4 . 34 Use of Acceleration Steps 

This method is a modification of steepest descent used to avoid 
the ravine problem. It uses steepest descent steps with occasional 
acceleration steps given by 

6 ^k = Pk C 5k ' *k-2> 

where = multiplier chosen for step-size control. 

These acceleration steps have the effect of moving along the axis of 
the zig-zag pattern shown in Figure 4.1 when the ravine problem is 
encountered. A modification of this procedure using alternating 
steepest descent and acceleration steps is given [34] . 
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4 ■ 4 Second-Order Techniques 


In finding a solution z=a to an equation f(z)=a, an iterative 
technique which acts such that 

z k+l ' a = ^V 01 ) 2 kk^ 1 (4.70) 

is known as a second-order technique. The technique used in this 
thesis is a generalized Newton-Rapheson iteration. It may be used 
only if the minimum value of the cost function J is zero since it 
searches for a zero of the u vector. 

Using the notation of (4.1) and (4.3), it is possible to con- 
struct the matrix equation 


-u = H 6s 


(4.71) 


where H is a matrix built up of the (3x3) sub -matrices H . according 

- 1 - » 3 

to the rules 


1) For i = odd. 


H i-l,i " D i + 1 


(=0 for i=l) 


H i,i " B i+1 ‘ P i+1 ‘ C i+1 ‘ D i+2 


H i+l,i " " A i+1 “ B i+2 ' P i+2 + C i+2 


H i+2,i ~ A i+2 


(=0 for i=2N-3) 


2) For i = even 


H i-l,i - - D i + 1 

H i,‘i = C i+1 ’ P i+1 * B i+1 


H i+l,i A i+1 


(=0 for i»2N-4) 


(4.72) 

(4.73) 

(4.74) 

(4.75) 

(4.76) 

(4.77) 

(4.78) 
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3) For all H. - not covered by (4 . 72) - (4 . 78) 
1 > 1 


0 0 



0 

0 


0 

0 

0 


(4.79) 


These matrices are identical to the bracketed terms in (4.38) and 

(4.40). The subscripts are offset by 1 since the terms in u and Ss^ 

run from 2 to 2N-3 while the subscripts of H run from 1 to 2N-4. 

1 > J 

Solving (4.71) yields the variation which would cancel the velocity 
mis-match vectoT if the linearization were an exact process 


6s = - H' 1 u 


(4.80) 


Since the linearization yields only a local approximation to the cost 
function surface, this process is iterative. However, second-order 
convergence is achieved. 

The matrix H to be inverted in (4.80) is of dimension 6(N-2) where 
N may be quite large. The inversion may be considerably simplified by 
noting that H is a banded matrix (i.e. has non-zero terms only on or 
near its main diagonal). Using this property solution (4.80) -may be 
obtained by solving the set of linear equations (4.71) using a Gauss- 
Jordan method which stores and manipulates only the non-zero elements. 
This allows the storage and computational time to be reduced consid- 
erably. The computer subroutine used is DGELB from [33]. 

4 . 5 Application of Inequality Constraints 

At a planet,- the coordinates of the entry and exit points are 
given by 
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(4.81) 


entry point = 



exit point = x k+1 


A W 

L At k + 1 


(4.82) 


Unit vectors centered m the- planet in the direction of the entry and 
exit points are given by 


'COS0 k 

cos4> k 

5in0 k 

COS(j»j. 

_simj> k 



ik+1 


cos(0 k + A0 k+1 )cos(^ lc + A4> k+1 ) 
S in(e k + A0 k+1 ) COS + AtJ> k+;l ) 

silica + A4> k+1 ) 


(4.83) 


(4.84) 


Assuming that these unit vectors approximate the asymptotes of the 
swmgby hyperbola, the turn angle v of the hyperbola is given by 


v = it - cos 


-1 


[ ik • ik+l ] 


(4.85) 


= TT 


cos 


[cos (Aej, + 1 )cos^> k cos C ( f > i c 4 -A4> k+1 ) 


+ sin^sinC^+A*^)] 


(4.86) 
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See Figure (4.2) for an illustration of this angle. The maximum 
turn angle possible for a minimum periapse radius of rp A is given by 


where 


= 2 sin 
M 


-1 


1 + r PA v I 


y i 


(4.87) 


vu = gravitational parameter of the planet 
Vj = hyperbolic excess velocity 


-H^k ^H,k ' 2 r 


s.i 


1 

1 


(4.88) 


r . = radius of planetary sphere of influence 
s , 1 


Before the acceptance of a step <Sx^ and | Sxj. + ^ for a planet's 
entry and exit points, the values of v and are calculated for the 
projected new values 


k = *k * % 
k* i = -k+a + 6 *k + i 


(4.89) 


If v>v M , the components of <Sxj, + 1 are set equal to zero and the 
components of 6x^ are left unchanged. This procedure is repeated for 
each swingby on every iteration. The criterion for the minimum periapse 
distance is 


r P,A ~ 


(4.90) 
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where 

r = equatorial radius of the planet 

C 

= constant multiplier (usually 1.1) 

4 . 6 Behavior of Minimization Techniques 

In general, the techniques discussed m sections 4.3 and 4.4 
behaved as follows; 

1) The steepest descent technique suffered badly from the 
ravine problem. It was not used extensively. 

2) The optimum gradient step (OGS) dealt with the ravine 
problem fairly well and also solved the step-size con- 
trol difficulties. It was adapted easily to the con- 
strained case. The OGS was the most useful first-order 
technique employed. 

3) The conjugate gradient method proved to be slightly 
faster than the OGS. It was not adaptable to the con- 
strained problem. Thus it was employed mainly as a check 
on the performance of the OGS. 

4) The acceleration steps showed no noticeable improvement 
over the OGS. They were not used extensively. 

5) The generalized Newton-Rapheson technique showed second- 
order convergence near the solution. However, it some- 
times diverged when started too far from the solution 
point. Also, it was applicable only when the minimum of 
the cost function was zero. 

In general, the first-order techniques showed a diminishing speed 
of convergence as they approach the solution point. They were useful 
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mainly to reach the neighborhood of the minimum. If applicable (i.e. 
if an admissible minimum equal to zero exists] the second-order tech- 
nique was used to converge to the solution. Otherwise the OGS technique 
was used for the whole process. 

4 . 7 Application of the Trajectory Segment Matching Technique 

The procedure for application of the trajectory segment matching 
procedure is as follows. 

1) Using the cost function given in (2.13) and (2.14) con- 
verge to a minimum of the cost for the advanced patched 
conic model. 

2) For this solution, compute the velocity offsets as 
described in Chapter 3. 

3) Using the cost function given in (3.31) and (3.32) con- 
verge to a minimum of the cost for the perturbed conic 
model. During this iteration, the offsets 6Vg ^ and 
8v^ ^ are considered to be constants and thus remain 
unchanged. 

4) For this new solution, re-compute the velocity offsets. 

If they have changed measurably, return to step 3. 
Otherwise, terminate the process. 

The result is a sequence of perturbed conic legs matched in position, 
velocity, and time at the sphere of influence entry and exit points. 

The accuracy of this sequence may be checked by determining a set of 
exact trajectory legs (corresponding to the same initial and final con- 
ditions on position and time) by numerical integration. In all cases 
examined, the perturbed conic legs have been close enough that the 
numerical trajectory integration has converged to the desired boundary 
conditions within two or three iterations. 

Examples of the use of this technique are given in the next three 
chapters . 
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-a 


Figure 4.1 Behavior of Steepest Descent Technique in a “Ravine 







Chapter 5 

Dual Planet Reconnaissance Trajectory Example 

5 . 0 Chapter Summary 

The dual planet reconnaissance trajectory is described. The coor- 
dinates of the sphere of influence entry and exit points and the orbit- 
al elements of the individual trajectory legs are given. The perfor- 
mance of the trajectory segment matching technique is discussed. The 
results of the perturbed conic analysis are tabulated and the analytic 
predictions are compared with numerically integrated trajectory legs. 
The resulting comparison shows that the trajectory may be predicted to 
within a total correction of 0.2263 m/sec. 

5 . 1 Description of the Mission 

This mission employs a free-fall trajectory (see Figure 5.1) which 
leaves Earth, makes a close pass first to Venus and then to Mars, and 
then returns to Earth. It would be possible for a reconnaissance mis- 
sion of both Mars and Venus by manned or unmanned spacecraft. In the 
unmanned case the return to earth would allow recovery of high res- 
olution photographs taken during the swingbys. This mission is de- 
scribed in [3] . A method for finding these dual planet swingbys using 
the simple patched conic model is given in [10] . 

The final trajectory is specified by the following set of plane- 
tary sphere of influence (SOI) entry and exit points. The SOI radii 


are given in 

Appendix C . 




Point 

Julian Date 



Location 

1 

2441478.80000 

Earth 

SOI 

exit point (launch) 

2 

2441634.11977 

Venus 

SOI 

entry point 

3 

2441637.99955 

Venus 

SOI 

exit point 

4 

2441787.28715 

Mars 

SOI 

entry point 

5 

2441792.32920 

Mars 

SOI 

exit point 

6 

2441949.20000 

Earth 

SOI 

entry point (arrival) 
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The spherical coordinates of the SOI entry and exit points given m a 
planetocentric ecliptic coordinate system (see Appendix A) are 


Point 

Azimuth 

Elevation 

1 

142.800° 

-19.600° 

2 

332.249° 

-2.637° 

3 

177.796° 

1.716° 

4 

43.066° 

4.118° 

5 

222.458° 

3.885° 

6 

8.300° 

5.700° 


The complete set of cartesian coordinates of the SOI entry and exit 
points in the heliocentric and planetocentric coordinate frames plus 
the heliocentric coordinates of the planet are given in Appendix E. 

The total time of flight is relatively short (about 1.3 years) 
and divided approximately evenly among the heliocentric legs . 


Leg 

Time of Flight (days) 

Central Angle Traversed 

1-2 

155.31977 

244.268° 

2-3 

3.87978 

193.491° 

3-4 

149.28760 

122.855° 

4-5 

5.04205 

184.410° 

5-6 

156.87030 

85.562° 


The legs (2-3) and (4-5) are planetocentric hyperbolas, while legs 
(1-2), (3-4) and (5-6) are heliocentric ellipses. The orbital elements 
of these legs are 


Leg 

Semi-Maj or 

Axis 

Eccentricity 

Inclination 

1-2 

0.80837 

a.u. 

0.25644 

3.348° 

2-3 

0.72733 

r P 

4.28637 

3.053° 

3-4 

1.07057 

a.u. 

0.37045 

3.290° 

4-5 

0.24513 

r P 

13.00436 

94.337° 

5-6 

1.06599 

a.u. 

0.37479 

1.310° 
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The semi-major axis is given in astronomical units (a.u.) for the 
heliocentric legs and in planetary radii for the planetocentric legs. 

The vis-viva energy at earth needed to launch on this trajectory 
is 


C 3 - 18 378 km 2 /sec 2 

while the re-entry velocity on return to earth is 
v R = 15,639.128 m/sec 

The periapse radius (in kilometers and planetary radii), the 
periapse velocity (m km/sec), and the turn angle (as defined m 
Section 4.5) for the planetocentric legs are 

Leg Periapse Radius Periapse Velocity Turn Angle 

2-3 14,461.578 (2.3903) 10.904491 13. 491° 

4-5 10,034.548 (2.9427) 7.737818 4.410° 

The time of periapse is the mid-point between the entry and exit times 

on the sphere of influence. 

All of the parameters listed in this section are for the two -body 
trajectory legs used as reference orbits for the final perturbed conic 
analysis . 

5 . 2 Performance of Trajectory Segment Matching Techniques 

The trajectory segment matching technique is repeated through 
several cycles. The first of these uses only the advanced conic model. 
The initial conditions are obtained from [10] where the trajectory was 
determined using the simple patched conic model. All successive cycles 
employ the perturbed conic model with the velocity offsets re-calculated 
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for each cycle. The column labeled method indicates whether the 


optimum gradient step (OGS) 

or the generalized 

Newton- Raphe son (GNR) 

procedure is used for that 

step. The cost function is defined by 

(2.13-.14) for the first cycle and by (3.31-.32) 

for successive cycles 

First Cycle: Iteration 

Method 

Cost Function (Rm/sec] 

Initial Conditions 

- 

3.87296 

1 

OGS 

2.81187 

2 

OGS 

0.35187 

3 

GNR 

0.29405 (10 -3 ) 

4 

GNR 

0.34828 (10" 9 ) 

5 

GNR 

0.36653 (10" 17 ) 

Second Cycle: Iteration 

Method 

Cost Function 

Initial Conditions 

- 

0.35299 (10" 3 ) 

1 

GNR 

Q. 14899 (10 -9 ) 

2 

GNR 

0.17063 (10' 17 ) 

Third Cycle: Iteration 

Method 

Cost Function 

Initial Conditions 

_ 

0.65610 (10" 8 ) 

1 

GNR 

0.10061 (10" 16 ) 

Fourth Cycle: Iteration 

Method 

Cost Function 

Initial Conditions 

_ 

0.31691 (10" 12 ) 


The change from OGS to GNR occurred when the cost function fell 

-14 

below 0.5. A cycle was terminated when the cost went below 1.0 (10 ) 

and the whole procedure was ended when the cost at the initial point 
of a new cycle was less than 1.0 (10 . The total running time of 

this process on the IBM 360/65 was 48 sec. 

5.3 Perturbed Conic Results 


The first step in the perturbed conic model was the calculation 
of the perturbations in the initial and final position and velocity 



for each leg due to the disturbing accelerations. The magnitudes of 
these perturbations (in 1cm and m/sec) for the last cycle are 


Leg 

<$r( v 

$v(t x ) - 

6r(t 2 ) 

«v(t 2 ) 

1-2 

35,161.42 

39.727 

9331.04 

25.252 

2-3 

1,482.11 

26.470' 

1335.23 

23.669 

3-4 

11,780.23 

26.099 

2069.76 

3.544 

4-5 

242.58 

3.347 

236.88 

3.248 

5-6 

4,769.78 

3.698 

6545.68 

15.886 


The next step is the calculation of the position and velocity 
perturbations at the mid-point that eliminate the position pertur- 
bations at the initial and final points . These mid-point perturbation 
magnitudes (in km and m/sec) for the last cycle are 


Leg 


5v(t M ) 

1-2 

9197.48 

4.882 

2-3 

324.23 

50.435 

3-4 

6059.22 

0.448 

4-5 

9.31 

1.330 

5-6 

3618.21 

0.723 

The last step is 

the calculation of the offsets of the initial and 

final velocities 

. These offset magnitudes (in m/sec) 

for the last 

cycle are 



Leg 

6v(t x ) 

«v(t 2 ) 

1-2 

34.286 

24.139 

2-3 

19.034 

14.558 

3-4 

22.322 

3.211 

4-5 

2.221 

2.183 

5-6 

2.997 

14.982 
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S . 4 Comparison with Numerically Integrated Results 


The accuracy of the analytic procedure is evaluated by comparison 
with numerically integrated trajectory legs determined using the 
following procedure: 

1) The initial conditions for the numerical integration are 
set equal to the conic position and the perturbed conic 
velocity at the initial time t^, 

“ IqC^) 

vCV = v Q (t ; ) + 6v(t 1 ) (5.1) 

[~ denotes values on the numerically integrated tra- 
jectory] 

,2) The N-body equations of motion and the differential 
equations for the N-body stqte transition matrix are 
integrated forward to the final time t 2 . The errors in 
position and velocity at the final time are given by 

Ar.(t 2 ) = £(t 2 ) ' £.q C^) (5.2) 

Av(t 2 ) = v(t 2 ) - {v 0 (t 2 ) + Sv(t 2 )] (5.3) 

The numerical integration routine is described in [26] . 

3) Using the N-body state transition matrix between t^ and 
t 2 , the change in v(t^) needed to eliminate Ar(t 2 ) is 
calculated. 




Steps 2) and 3) are repeated until Ar(t 2 ) is driven to zero. If 
the initial guess is close this happens quite rapidly. For this example, 
A£.(t 2 ) was reduced to under 10 km in two iterations. 

The quantities chosen to indicate the accuracy of the analytic 
techniques are the magnitude of the velocity error at t-^ 

A lCt 1 ) = vft^ - [v 0 (t 1 ) + dvCtjM, (S • 4] 

the magnitude of the velocity error at t^ 

A lCt 2 ) = v(t 2 ) - [v Q Ct 2 ) + 6v(t 2 D], (5-5) 


and the position error at due to the use of .the analytically deter- 
mined velocity at for the numerical integration 

Ar(t 2 ) = " 10^2^ ( S ' 6 J 

for 

i(tl) = + <Sv(t 1 )' . 


These values 

(in km 

and m/sec) are 


Leg 

Av(t 1 ) 

av(t 2 ) 

Ar(t 2 ) 

1-2 

.1250 

.0450 

5327.23 

2-3 

.0030 

.0091 

9.71 

3-4 

. 0430 

.0041' 

1306.60 

4-5 

. 0028 

.0030 

1.62 

5-6 

.0017 

.0173 

15.39 

In order to 

fly the 

trajectory, it would be 

necessary to apply 


velocity corrections at each of the SOI entry and exit points to 
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eliminate these errors. The magnitudes of these corrections (in m/sec) 


are 


Point 


Correction 


1 

2 

3 

4 

5 

6 


0.1250 

0.0448 

0.0340 

0.0022 

0.0030 

0.0173 


The total correction needed is 


Av = 0 . 2263 m/sec. 


It is also possible to compare the periapse conditions of the 
numerically integrated trajectory with those of the two-body reference 
trajectory used m the perturbed conic analysis. The differences are 


Leg 

Ar 

TT 

Av 

IT 

At: 

IT 

2-3 

8.3576 km 

8.0255 m/sec 

-29.142 sec 

4-5 

0.6106 km 

1.1715 m/sec 

-1.223 sec 


where 

t = time of periapse for two-body reference orbit 
0 

t^ = time of periapse for numerically integrated trajectory 


ir = 

IT — TTq^TTq^ 
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5.5 Discussion of Results 


The analytic technique has provided an accurate description of the 
legs of the dual planet reconnaissance trajectory. The velocity errors 
described in the last section could be eliminated by running the tra- 
jectory segment matching procedure for a few more cycles using offsets 
calculated by numerical integration rather than by the analytic tech- 
niques of Chapter 3. Instead, it is probably easier to absorb these 
errors in the mid-course corrections made on approach and departure 
from each planet. 
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Figure 5.1 Dual Planet Reconnaisanpe Trajectory 
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Chapter 6 

Grand Tour Trajectory Example 

6 . 0 Chapter Summary 

The Grand Tour mission is described. The coordinates of the 
sphere of influence entry and exit points and the orbital elements 
of the individual trajectory legs are given. The performance of the 
trajectory segment matching procedure is discussed. The results of 
the perturbed conic analysis are tabulated and the analytic predic- 
tions are compared with numerically integrated trajectory legs. This 
comparison shows that the trajectory may be predicted to within a 
total correction of 2.652 m/sec. The accuracy of the model for the 
disturbing acceleration during the planetocentric phase is considered. 

6 ■ 1 Description of the Mission 

This mission employs a free-fall trajectory (see Figures 6.1 
and 6.2) which leaves Earth and makes successive close passes to 
Jupiter, Saturn, Uranus, and Neptune. The configuration of the planets 
necessary for the Grand Tour occur only once every 179 years with the 
next opportunity occurring m the period 1975-1981. Descriptions of 
this type of mission may be found m [16], [17], and [20] where the 
simple patched conic model is used to determine launch windows and 
approximate trajectory parameters. The specific trajectory chosen for 
this chapter leaves Earth in October, 1978. 

The final trajectory is specified by the following set of plane- 
tary sphere of influence (SOI) entry and exit points. The SOI radii 


are listed 

in Appendix C. 


Point 

Julian Date 

Location 

1 

2443787.00000 

Earth SOI exit point (launch) 

2 

2444291.61927 

Jupiter SOI entry point 
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Point 

Julian Date 


Location 

3 

2444457.85596 


Jupiter SOI exit point 

4 

2444930.62043 


Saturn SOI entry point 

5 

2445126.30248 


Saturn SOI exit point 

6 

2446481.99628 


Uranus SOI entry point 

7 

2446638.44605 


Uranus SOI exit point 

8 

2447720.00000 


Neptune SOI entry point 




(arrival) 

The launch date is 

October 6, 1978 

while 

the arrival date is July 13, 

1989 (which luckily 

occurs on a Thursday 

in that year) . 

The spherical 

coordinates of 

the SOI 

entry and exit points given 

in a planetocentric 

ecliptic coordinate system (see Appendix A) are 


Point 

Azimuth 

Elevation 

1 

101.424° 

8.720° 

2 

312.920° 

1.093° 

3 

180.066° 

4.253° 

4 

16.039° 

-3.433° 

5 

281.565° 

-2.502° 

6 

94.166° 

1.665° 

7 

291.444° 

2.975° 

8 

114.270° 

-2.812° 


The complete set of cartesian coordinates of the SOI entry CtliU CAi L. 
points in the heliocentric and planetocentric coordinate frames plus 
the heliocentric coordinates of the planets at encounter are given in 
Appendix B. 

The total flight time for the trajectory is short (about 10.77 
years) considering the distance covered. The time of flight and central 
angle traversed for the individual legs are 
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Time of Flight (days) Central Angle Traversed 


1-2 

504.61927 

144.498° 

2-3 

166.23669 

205.609° 

3-4 

472.76447 

25.902° 

4-5 

195.68205 

223.050° 

5-6 

1355 . 69380 

54.438° 

6-7 

156.44977 

189.035° 

CO 

1 

1081. 5539S 

17.392° 


The legs (2-3), (4-5), and (6-7) are planetocentric hyperbolas 
about Jupiter, Saturn, and Uranus respectively. The leg (1-2) is a 
heliocentric ellipse while legs (3-4), (5-6) and (7-8) are heliocentric 
hyperbolas. Each swingby adds energy to the trajectory, keeping it 
above solar escape energy after the Jupiter swingby. The orbital ele- 
ments for the trajectory legs are 


Leg 

Semi-Major Axis 

Eccentricity 

Inclination 

1-2 

4.61480 a.u. 

0.78328 

2.383° 

2-3 

17.22512 r p 

2.31355 

6.881° 

3-4 

28.43593 a.u. 

1.14287 

2.857° 

4-5 

5.23592 r p 

1.46490 

4.412° 

5-6 

3.77866 a.u. 

3.54820 

2.836° 

6-7 

1.10298 r p 

6 .,36816 

15.110° 

7-8 

3.06068 a.u. 

5.41387 

2.821° 


The semi-major axis is given m astronomical units (a.u.) for the 
heliocentric legs and in planetary radii for the planetocentric legs. 

The vis-viva energy at earth needed to launch on this trajectory 
is 

C 3 = 101.520 km 2 /sec 2 
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The periapse radius (in kilometers and planetary radii), the 
periapse velocity (in km/sec) and the turn angle (as defined in Section 
4.5) for the planetocentric legs are 


Leg 

Periapse Radius 

Periapse Velocity 

Turn Angle 

2-3 

1,615,506.286 (22.62614) 

16.121540 

25.609° 

4-S 

147,025.929 (2.4342) 

25.207545 

43.050° 

6-7 

139,143.535 (5.9210) 

17.530135 

9.035° 


The time of periapse is the mid-point between the entry and exit times 
on the SOI. 

All the parameters listed in this section are for the two-body 
trajectory legs used as reference orbits for the final perturbed conic 
analysis . 

6 , 2 Performance of Trajectory Segment Matching Technique 

The traj’ectory segment matching procedure is repeated through 
several cycles. The first of these uses the advanced patched conic 
model. The initial conditions are obtained from [16] where the tra- 
jectory was determined using the simple patched conic model. All suc- 
cessive cycles employ the perturbed conic model with offsets re-calcu- 
lated for each cycle. The column labelled method indicates whether the 
optimum gradient step (OGS) or the generalised Newton-Rapheson (GNR) 
procedure is used for that step. The cost function is defined by 
(2.13-.14) for the first cycle and by (3.31-.32) for successive cycles. 

2 ? 

First Cycle: Iteration Method Cost Function (krtT/sec ) 


Initial Conditions 

- 

4.04945 

1 

OGS 

3.96201 

2 

OGS 

3.89372 

3 

OGS 

3.79476 

4 

OGS 

3.76557 
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First Cycle: 

Iteration 

Method 

2 2 

Cost Function (km /sec ) 

5 


OGS 

3.72571 

92 


OGS 

1.47623 

93 


OGS 

1.47268 

94 


OGS 

1.46821 

95 


GNR 

0.47856 (10' 2 ) 

96 


GNR 

0.23789 (10 -7 ) 

97 


GNR 

0.80956 (10" 13 ) 

98 


GNR 

0.10929 (10" 18 ) 

Second Cycle 

• Iteration 

Method 

Cost Function 

Initial 

Conditions 

- 

0.10898 (10' 1 ) 


1 

GNR 

0.10408 (10 -5 ) 


2 

GNR 

0.13908 (lO" 11 ) 


3 

GNR 

0.32573 (10~ 17 ) 

Third Cycle: 

Iteration 

Method 

Cost Function 

Initial 

Conditions 


0.10560 (10 -5 ) 


1 

GNR 

0.10959 (10 -11 ) 


2 

GNR 

0.19818 (10' 17 ) 

Fourth Cycle 

. Iteration 

Method 

Cost Function 

Initial 

Conditions 

_ 

0.15182 (10‘ 9 ) 


1 

GNR 

0.24154 (10~ 15 ) 

Fifth Cycle: 

Iteration 

Method 

Cost Function 

Initial 

Conditions 

. 

0.15638 (10' 11 ) 


The change from OGS to GNR occurred when the cost function fell 
below 1.47. A cycle was terminated when the cost fell below 1.0 (lO -1 ^) 
and the whole procedure was ended when the cost at the beginning of 
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a new cycle was less than 1.0 (10~ 10 ). The total running time on the 
IBM 360/65 was about 172 seconds. 

6 . 5 Perturbed Conic Results 

The first step in the perturbed conic model was the calculation 
of the perturbations m the initial and final position and velocity 
for each leg due to the disturbing accelerations. The magnitudes of 
these perturbations (in km. and m/sec) for the last cycle are 


Leg 

6r{t 1 ) 

SvCtp 


6r(t 2 } 

6v(t 2 ) 

1-2 

144,135.26 

20.974 


680,081.06 

115.831 

2-3 

350,124.98 

146.486 


317,175.98 

126.799 

3-4 

681,236.05 

117.365 


90,828.50 

7.226 

4-5 

106,229.23 

38.040 


51,438.66 

18.086 

5-6 

847,681.45 

46.614 


464,247.07 

12.626 

6-7 

8,764.69 

3.921 


7,684.89 

3.384 

7-8 

246,233.21 

11.917 


222,593.00 

10.058 

The 

next step is the 

: calculation 

of the position and velocity 

perturbations at the mid- 

point that 

eliminate the 

position perturbations 

at the initial and final 

points. The magnitudes of 

these mid-point per- 

turbations (in km and m/sec) for the 

last cycle are 


Leg 


^(t M ) 



dv(t M ) 

1-2 


276,302.85 



13.418 

2-3 


14,105.26 



56.324 

3-4 


367,398.01 



14.229 

4-5 


174,611.87 


12 

,147.092 

5-6 


616,749.98 



3.755 

6-7 


1,673.59 



29-264 

7-8 


215,439.75 



2.071 

The 

last step is the 

calculation 

of the offsets of the initial 
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and final velocities. These offset magnitudes (in m /sec) for the last 


cycle are 


Leg 

SvCt x ) 

6v(t 2 ] 

i 

1 

16.800 

92.557 

2-3 

94.702 

81.347 

3-4 

100.210 

10.165 

4-5 

34.432 

20.439 

5-6 

42.388 

14.480 

6-7 

2.820 

2.104 

7-8 

10.560 

9.195 


6 , 4 Comparison with Numerically Integrated Results 

The accuracy of the analytic procedure is evaluated by comparison 
with numerically integrated trajectory legs using the procedure de- 
scribed in Section 5.4. The quantities used to evaluate the accuracy 
of the analytic technique are the errors in the velocity offsets at 
the initial and final times and the position error at the final time 
due to the use of the analytically determined initial velocity. The 
magnitudes of these quantities (in km and m/sec) are 


Leg 

Avftj) 

Av(t 2 ) 

Ar,(t 2 ) 

1-2 

0.305 

0.139 

8185.05 

2-3 

0.622 

0.438 

8955.75 

3-4 

0.353 

0.089 

15,348.52 

4-5 

0.335 

0.184 

239,662.60 

5-6 

0.052 

0.095 

_ 6444.35 

6-7 

0.0081 

0.0071 

1904.25 

7-8 

0.095 

0.095 

8833.50 

In order to 

fly the trajectory, 

it would be 

necessary to apply 


velocity corrections at each SOI entry and exit point to eliminate the 
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above errors. The magnitudes of these corrections (in m/sec) are 
Point Correction 

1 0.30S 

2 0.667 

3 0.787 

4 0.424 

5 0.179 

6 0.103 

7 0.092 

8 0.095 

The total correction needed is 
Av = 2.652 m/sec. 

The differences in the periapse conditions between the numerically 
integrated trajectory legs and the two-body reference values for the 
planetocentnc legs are 


Leg 

Ar 

TT 


AV^ 

TT 

At 

TT 

2-3 

1025.172 

km 

29.075 m/sec 

18.5368 minutes 

4-5 

40.068 

km 

5.083 m/sec 

-116.5053 minutes 

6-7 

8.162 

km 

1.216 m/sec 

-1.5907 minutes 


6 . 5 Discussion of Results 

The assumption that the disturbing acceleration due to other 
planets has neglible effect during the planetocentric legs of the tra- 
jectory does not prove to be as good for this example as it did for the 
dual planet reconnaissance trajectory. This shows up most strongly for 
the planetocentric leg about Saturn due to the long time inside the 
sphere of influence (196 days) and the proximity of Jupiter (about 
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5 a.u. away). To evaluate the effects of the other planets, the 
numerical integration of the Saturn planetocentric leg was repeated 
using the sun as the only disturbing body. The results are 


Av(t^) = 0.026 m/sec (0.335 m /sec) 

Av(t 2 ) = 0.004 m/sec (0.184 m/sec) 

Ar(t 2 ) = 1704.871 km (239,662.60 km) 

Ar^ = 37.827 km (40.068 km) 

Av^ = 4.918 m/sec (5.083 m/sec) 

At^ = -115.4210 minutes (-116.5053 minutes) 


The corresponding figures for the numerical integration of the same 
leg using all the other planets as well as the sun as disturbing 
bodies are shown in parenthesis. 

The large shift m periapse time at Saturn shows up in large 
values of 6r(t^) and 6v(t^) given in Section 6.3. Position and veloc- 
ity change quite rapidly near periapse for a hyperbola. 

The general accuracy of the analytic procedure as applied to 
this trajectory is quite adequate. The total correction needed (2.652 
m/sec) may be quite easily absorbed into the mid-course correction 
allowance. Based on the results described above, about 401 of the 
total correction magnitude could be eliminated by using a perturbed 
conic model which includes the disturbing accelerations of the other 
planets during the planetocentric phases . 
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Figure 6-J Grand Tour Trajectory 


98 






Figure 6.2 Grand Tour Trajectory Segment 
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Chapter 7 

Periodic Trajectory Example 

7 . 0 Chapter Summary 

The periodic trajectory example is described. The coordinates 
of the sphere of influence entry and exit points and the orbital 
elements of the individual trajectory legs are given. The performance 
of the trajectory segment matching procedure is discussed. The re- 
sults of the perturbed conic analysis are tabulated and the analytic 
predictions are compared with numerically integrated trajectory legs. 
The solution found is not a free-fall trajectory but requires a total 
impulse of 220.534 m/sec applied at the various sphere of influence 
(SOI) entTy and exit points . The total error in the calculation of 
the trajectory legs amounted to 38.950 m/sec which would also be 
applied at the SOI entry and exit points. The sources of these errors 
are discussed. 

7 . 1 Description of Mission 

The term periodic orbit is used here to mean an interplanetary 
free-fall trajectory which shuttles bach and forth between two planets. 
Once this orbit is established, it continues indefinitely with only 
minor guidance corrections. The existence of such trajectories was 
first explored by Hollister in [22]. Using the simple patched conic 
model, a general search procedure for periodic trajectories was de- 
veloped and three types of trajectories connecting Earth and Venus 
were discovered. This work was extended by Menning [23] who found 
additional Earth-Venus trajectories and by Rail [25] who examined 
Earth-Mars periodic trajectories. The guidance requirements for the 
Earth-Venus trajectories were examined by Hickman in [24]. 

The trajectory chosen for closer examination using the technique 
developed m this thesis is a segment of Periodic Orbit I given in 
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[22]. It consists of the sequence 

1) Earth to Venus transfer 

2) Direct Return to Venus 

3) Direct Return to Venus 

4) Venus to Earth transfer 

5) Direct Return to Earth 

6) Earth to Venus transfer 

The direct return legs are trajectories which return to the launch 
planet one planetary year later. This sequence (illustrated in Figures 
7.1 to 7.3) forms a part of a larger (16 year) orbit which repeats 
indefinitely between Earth and Venus. 

The final trajectory is specified by the following set of plane- 


tary sphere of influence (SOI) entry and exit points. The SOI radii 


are given m 

Appendix C. 




Point 

Julian Date 



Location 

1 

2440811. 2980S 

Earth 

SOI 

exit point (launch) 


• 

. - <* 

i j 

V - i ‘ 

2 

2440968.17351 

Venus 

SOI 

entry point 

3 

i 1. ^ 

2440974.71063 

Venus 

SOI 

exit point 

1 IlJ ^ 

4 

2441192 .53554' 

f j *"j 

i * ' ’ 

Venus 

SOI 

entry point 

V ' 

2441199.05152 

I ^ ^ ^ 

Venus 

SOI 

exit point 

6 ' 

2441416.81361’ 

r - » 

Venus 

■* # . 

SOI 

. it * i 

entry point 

" f , 

■ i 1 

•* H ' ! 

t ■* ■* 

- - * 1 

7 

2441423.35856 

Venus 

SOI 

exit point 

l 

7 ' a 

*, « • 

, cr 


8 

2441587.41963 

Earth 

SOI 

entry point 

,.U)J 

-u-' -,y v.f ... 

~ i 

! 


9 

2441598.42073 

Earth 

SOI 

exit point 

"* , **J j 

Vo 

1 _• J f ^.t ; > 

2441951.97897 

Earth 

SOI 

entry point 

*- i 

11* 

2441962.96593 

Earth 

SOI 

} i L. 

exit point 

12 

2442123.36905 

Venus 

SOI 

entry point (arrival) 


The launch date is on August 11, 1970 with the arrival date on March 
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13, 1974 (which is a Wednesday). 

The spherical coordinates of the SOI entry and exit points in a 


planetocentric 

ecliptic coordinate system 

(see Appendix A) : 

Point 

Azimuth 

Elevation 

1 ' 

-68.845° 

-71.474° 

2 

35.165° 

-52.462° 

3 

157.616° 

18.414° 

4 

336.896° 

-16.982° 

5 

154.209° 

-54.071° 

6 

334.844° 

55.489° 

7 

90.723° 

-121.732° 

8 

76.408° 

-48.250° 

9 

3.803° 

41.428° 

10 

183.950° 

-40.474° 

11 

314.174° 

-42.542° 

12 

151.149° 

-65.069° 


The complete set of cartesian coordinates of the SOI entry and exit 
points in the heliocentric and planetocentric coordinate frames plus 
the heliocentric coordinates of the planets at encounter are given in 
Appendix E. 

The total flight time for the trajectory is 3.59 years which is 
about one-fifth of the basic periodic orbit period. The time of flight 
and central angle traversed for the individual legs are 


heg 

Time of Flight (days) 

Central Angle Traversed 

1-2 

156.87546 

196.498° 

2-3 

6.53712 

208.873° 

3-4 

217.82491 

348.950° 

4-5 

6.51598 

216.238° 
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Leg 

Time of Flight (days) 

Central Angle Traversed 

5-6 

217.76209 

348.831° 

6-7 

6.54495 

208.537° 

7-8 

164.06107 

198.183° 

8-9 

11.00109 

215.694° 

9-10 

353.S5824 

348.607° 

10-11 

10.98699 

227.708° 

11-12 

159.40312 

192.359° 

The legs 

(2-3), (4-5), (6-7), (8-9) 

and (10-11) are planetocentric 


hyperbolas while the legs (1-2), (3-4), (5-6), (7-8), (9-10), and 
(11-12) are heliocentric ellipses. The orbital elements of these legs 
are 


Leg 

Semi-Major Axis 

Eccentricity 

Inclination 

1-2 

0.86448 

0.17688 

7.058° 

2-3 

2.15499 

2.07092 

126.095° 

3-4 

0.72333 

0.13795 

4.093° 

4-5 

2.14851 

1.69163 

91.594° 

5-6 

0.72333 

0.07569 

7.270° 

6-7 

2.16000 

2.09324 

70.943° 

7-8 

0.86362 

0.16973 

6.908° 

8-9 

3.26454 

1.71393 

59.504° 

9-10 

0.99995 

0.09105 

5.438° 

10-11 

3.27180 

1.35186 

115.417° 

11-12 

0.86564 

0.16374 

6.373° 

semi-maj or 

axis is given m 

astronomical units 

for the heliocentrii 


legs and in planetary radii for the planetocentric legs. 

The vis-viva energy at earth needed to launch on this trajectory 
is 

C 3 = 18.427 hm 2 /sec 2 
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The periapse radius fin kilometers and planetary radii) , the 
periapse velocity (in km/sec) and the turn angle (as defined m Sec- 
tion 4.5) for the planetocentric legs are 


Leg 

Periapse 

Radius 

Periapse Velocity 

Turn Angle 

2-3 

13,962.344 

(2.30783) 

8.45835 

28.873° 

4-5 

8,990.130 

(1.48597) 

9.86860 

36.238° 

6-7 

14,286.449 

(2.36140) 

8.39219 

28.537° 

8-9 

14,865.250 

(2.33065) 

8.52453 

35.694° 

10-11 

7,342.533 

(1.15120) 

11.29118 

47.708° 

time of 

periapse is the 

mid-point 

between the entry and 

exit points 


on the sphere of influence. 

All of the parameters listed in this section are for the two-body 
trajectory legs used as reference orbits for the final perturbed conic 
analysis . 

7 . 2 Performance of Trajectory Segment Matching Technique 

No free-fall solution was found for the trajectory sequence chosen 
for detailed examination. Since no zero of the cost function existed, 
the generalized Newton-Rapheson (GNR) technique could not be used. All 
attempts at application of GNR resulted in a rapid divergence from the 
minimum. Thus, the trajectory segment matching procedure was restricted 
to the first-order optimum gradient step (OGS) technique. This tech- 
nique was repeated through five cycles. The first of these used the 
advanced patched conic model alone with initial conditions obtained 
from [22] and [24]. All successive cycles employed the perturbed conic 
model with the velocity offsets re-computed for each cycle. The cost 
function is defined by (2.13-.14) for the first cycle and by (3.31-.32) 
for successive cycles. 
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Method 


2 2 

Cost Function (1cm /sec ) 


First Cycle : Iteration 


Initial Conditions 

- 

0.61931 

1 

OGS 

0.41670 

2 

OGS 

0.25163 

3 



198 

OGS 

0.0006968 

199 

OGS 

0.0006964 

200 

OGS 

0.0006961 

Second Cycle: Iteration 

Method 

Cost Function 

Initial Conditions 

- 

0.023412 

1 

OGS 

0.013790 

2 

OGS 

0.013629 

3 

OGS 

0.013415 

49 

OGS 

0.012933 

SO 

OGS 

0.012926 

Third Cycle: Iteration 

Method 

Cost Function 

Initial Conditions 

- 

0.012836 

1 

OGS 

0.012824 

2 

OGS 

0.012816 

79 

OGS 

0.012250 

80 

OGS 

0.012241 

Fourth Cycle: Iteration 

Method 

Cost Function 

Initial Conditions 

- 

0.012155 

1 

OGS 

0.012140 

2 

OGS 

0.012131 
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Fourth Cycle: Iteration 

Method 

Cost Function 

99 

OGS 

0.011551 

100 

OGS 

0.011545 

Fifth Cycle: Iteration 

Method 

Cost Function 

Initial Conditions 


0.011491 


As the iteration approached a minimum the rate of convergence 
slowed down. The procedure was terminated at the fifth cycle due to 

1) The small change m the solution point during the third 
and fourth cycles (the coordinates varied by 10 5 days 
and 10 3 degrees) . 

2) The small change in the cost function for each iteration 
(about 10" S to 10~ 6 per iteration in the fourth cycle). 

3) The small change in the cost function between the fourth 
and fifth cycles. 

6J = 5.4 (10' 5 ) 

4) Fluctuations m the components of the OGS weighting 
matrix indicating that the path of the iteration changes 
direction repeatedly. 

Since the trajectory sequence found is not a free-fall solution, 
a velocity impulse is required at each SOI entry and exit point. The 
magnitudes of these impulses (in m/sec) for the last cycle are 

Point Impulse Required 

1 

2 0.181 

3 9.935 
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Point 


Impulse Required 


4 

5 

6 

7 

8 
9 

10 

11 

12 


10.380 

69.168 

70.S61 

0.607 

2.240 

27.545 

27.402 

2.515 


Total Impulse 220.534 m/sec 

Since they are boundary points and need not be matched with any other 
trajectory segments, points 1 and 12 do not have impulses required. 


Total running time on the IBM 360/65 was about 450 seconds. 


7 . 3 Perturbed Conic Results 

The first step in the perturbed conic model is the calculation 
of the perturbations in the initial and final position and velocity 
for each leg due to the disturbing accelerations. The magnitudes of 
these perturbations (in km -and m/sec) for the last cycle are 


Leg 

6r(t 1 ) 

8v(t 1 ) 

Sr(t 2 ) 

Sv(t 2 ) 

1-2 

41,314.08 

40.979 

32,426.32 

40.633 

2-3 

2614.02 

26.988 

4091.55 • 

42.725 

3-4 

28,220.08 

45.670 

29,609.92 

46.557 

4-5 

4237.34 

44.509 

2906.04 

30.227 

5-6 

46,910.76 

49.855 

48,591.14 

50.337 

6-7 

2821.81 

29.843 

2253.46 

23.601 

CO 

1 

31,659.97 

40.791 

40,915.69 

38.892 

8-9 

3542.84 

21.977 

5319.98 

33.174 
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Leg 

SrC^) Svtt.^ 

«r(t 2 ) 

Sv(t 2 ) 

9-10 

48,817.55 43.424 

59,081.77 

46.252 

10-11 

5558.97 34.485 

4081.22 

25.052 

11-12 

40,283.06 39.210 

23,510.72 

35.092 

The next 

step is the calculation 

of the position and velocity 

perturbations 

at the mid-point that eliminate the 

position perturba- 

tions at the initial and final points. 

The magnitudes of these mid- 

point perturbations (m km and m/sec) 

for the last cycle aTe 

Leg 



Sv(t M ) 

1-2 

200,104.06 


13.023 

2-3 

4079.14 


810.701 

3-4 

1,096,815.12 


249.635 

4-5 

4147.93 


1694.559 

5-6 

1,415,883‘.04 


339.272 

6-7 

1170.26 


224.513 

7-8 

184,637.38 


10.833 

8-9 

4769.56 


1010.601 

9-10 

1,845,405.84 


271.268 

10-11 

5340.29 


3490.696 

11-12 

242,051.23 


12.239 

The last 

step is the calculation 

of the offsets of the initial 

and final velocities. These offset magnitudes (in m/sec) for the last 

cycle are 




Leg 

«V(tj) 


6v (t 2 ) 

1-2 

18.892 


47.168 

2-3 

23.658 


36.195 

3-4 

136.154 


137.953 

4-5 

36.352 


22.588 
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Leg 

5v(t 1 j 

6v(t 2 ) 

5-6 

179.485 

180.798 

6-7 

19.858 

17.587 

7-8 

36.508 

26.500 

8-9 

18.928 

25.041 

9-10 

168.229 

145.449 

10-11 

25.288 

19.274 

11-12 

37.285 

43.667 


7 . 4 Comparison with Numerically Integrated Results 

The accuracy of the analytic procedure is evaluated by comparison 
with numerically integrated trajectory legs using the procedure de- 
scribed in Section 5.4. The quantities used to evaluate the accuracy 
of the analytic technique are the errors in the velocity offsets at 
the initial and final times and the position error at the final time 
due to the use of the analytically determined initial velocity. The 
magnitudes of these quantities (in km and m/sec) are 


Leg 

Avftj) 

Av(t 2 ) 

Ar(t 2 ) 

1-2 

0.265 

0.397 

6474.90 

2-3 

0.019 

0.078 

307.75 

3-4 

4.091 

4.148 

3645.73 

4-5 

0.072 

0.020 

453.86 

5-6 

8.757 

8.830 

7159.38 

6-7 

0.030 

0.010 

360.25 

7-8 

0.279 

0.238 

7305.56 

8-9 

0.007 

0.052 

143.39 

9-10 

5.489 

5.574 

5829.96 

10-11 

0.045 

0.012 

1849.09 

11-12 

0.388 

0.240 

8008.43 
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In order to fly the trajectory, it would be necessary to apply 
additional velocity corrections at each SOI entry and exit point to 
eliminate the above errors. The following table lists the magnitudes 
(in m/sec) of the impulse required at each SOI entry and exit point 
as predicted by the analytic technique and the magnitudes (in m/sec) 
of the errors in these predictions as- determined by the numerical 


integrations . 

Point Impulse Required Error 

1 - 0.265 

2 0.181 0.410 

3 9.935 4.166 

4 10.380 4.218 

5 69.168 8.749 

6 70.561 8.856 

7 0.607 0.277 

8 2.240 0.239 

9 27.545 5.532 

10 27.402 5.612 

11 2.515 0.386 

12 - 0.240 

Totals 220.534 38.950 


The differences m-the periapse conditions between the numerically 
integrated trajectory legs and the two-body reference values for the 
planetocentric legs are 


Leg 

Ar 

IT 


Av,, 

7f 

At 

7T 

2-3 

18.01 

km 

11.351 m/sec 

7.997 minutes 

4-5 

19.81 

km 

6.743 m/sec 

-6.991 minutes 

6-7 

38.27 

km 

5.371 m/sec 

-2.307 minutes 

8-9 

66.57 

km 

14.038 m/sec 

9.291 minutes 

10-11 

38.36 

km 

31.352 m/sec 

7.879 minutes 
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7 . 5 Discussion of Results 


Singe the solution found is not a zero of the cost function, it 
may not be assumed to be a global minimum. The divergence of the GNR 
iteration indicates strongly that no zero of the cost function exists 
locally but this has not been proven. Similarly, the existence of 
other lower but non-zero minima has not been disproven. Some experi- 
mentation using different sets of initial conditions was conducted 
but no firm conclusions on the existence of multiple solutions were 
reached. 

The calculation of the direct return trajectories in this example 
proved to be a most severe test of the trajectory determination tech- 
nique. The initial and final velocity offsets required for these legs 
were an order of magnitude larger than those for previous trajectories. 
The source of these large offsets may be seen by examining the tra- 
jectories in Figures 7.1 to 7.3. The direct return trajectories ob- 
viously spend a large time in the vicinity of the planet with which 
they are associated leading to their somewhat extreme behavior. 

The errors in the offset calculations associated with the direct 
return trajectories are also about an order of magnitude larger than 
those for previous heliocentric trajectories. This is due to two effects 

1) The larger offsets put a greater strain on the linearity 
assumptions of the perturbed conic model. 

2) The nearness of the associated planet over a large part 
of the trajectory causes significant differences between 
the two -body state transition matrix used and the actual 
many-body state transition matrix. 

Of these two sources, the second is considered more significant. One 
point which is interesting to note is the fact that the final position 
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error due to using the analytically determined initial velocity is 
not significantly larger for the direct return trajectories. This 
indicates that the velocity errors for the direct returns are not 
m a critical direction. 

The general accuracy of the trajectory determination procedure 
as applied to the periodic trajectory example appears to be compatible 
with the results of the preceding examples. For applications for which 
the accuracy of the direct return leg calculations proved unacceptable, 
a third stage using numerical integration to determine the velocity 
offsets for these legs could be added to the solution procedure. 
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1 a. u . 

‘■ ' i - J 


Figure 7,1 Period ic Trajectory S 
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Figure 7.2 Periodic Trajectory Segment 
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1 a . u . 



Figure 7.3 Periodic Trajectory Segment 
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Chapter 8 

Summary, Conclusions, and Recommendations 
8 . 0 Summary and. Conclusions 

The trajectory targeting technique developed in this thesis is 
intended for the pre-mission calculation of reference trajectories 
for multiple swingby interplanetary trajectories . The primary objective 
has been to develop a technique which 

i) has a wide enough region of convergence such that 
the initial conditions for the trajectory determina- 
tion process may be derived from a simple patched 
conic mission analysis, 

ii) is largely analytic to minimize computational time 
required, and 

ni) is accurate enough to eliminate or at least signif- 
icantly reduce the need for numerical integration of 
traj ectories . 

The trajectory targeting technique developed is applied in the 
following manner. A simple patched conic model [consisting of a sequence 
of heliocentric conic arcs running from the center of one planet to the 
next matched in relative hyperbolic excess velocity at each planetary 
encounter) is used to calculate a set of initial conditions for an 
advanced patched conic model. This advanced patched conic model con- 
sists of a set of alternating heliocentric and planetocentric conic 
arcs joined at the planetary sphere of influence (SOI) . These arcs are 
specified by the position and time of the entry and exit points of the 
trajectory through the SOI of each planet encountered (i.e. the helio- 
centric arcs run from the exit point on one SOI to the entry point on 
the next SOI while the planetocentric arcs run from the entry point to 
the exit point of a single SOI) . Since the end points and the time of 


117 


flight for each arc are specified, the conic initial and final veloc- 
ities may be calculated by solving Lambert's Problem. 

Since the initial conditions for the SOI entry and exit points 
were determined using an approximate model, the conic arcs in the 
advanced patched conic model are not likely to match dynamically. 
Instead, velocity discontinuities occur at each SOI entry and exit 
point. Using the sum of the squares of the magnitudes of these veloc- 
ity mis-matches as a cost function, the next step is to formulate the 
problem of varying the entry and exit points and times to minimize 
the total velocity mis -match as a parameter optimization problem. The 
expression for the gradient of the cost function with respect to the 
entry and exit point coordinates may be determined analytically from 
the relations for the Lambert Problem. Then, by applying first- or 
second-order iteration techniques, the velocity mis -match may be 
reduced to a minimum, which will be zero for a free-fall trajectory. 

Once the velocity mis-match has been minimized for the advanced 
patched conic model the next step is to repeat the process using the 
perturbed conic model. Using the two-body conic arcs calculated for the 
advanced patched conic model as reference trajectories, perturbed conic 
segments (which include perturbations caused by the disturbing accel- 
erations of the planets on the heliocentric legs and the disturbing 
acceleration of the sun on the planetocentric legs) which pass through 
the same end points and times as the advanced patched conic segments 
are calculated analytically. These perturbed conic segments differ from 
the advanced patched conic segments by velocity offsets at the initial 
and final times. The cost function is now modified to include these 
offsets and the iteration procedure to minimize the velocity mis-match 
(now including the offsets) is repeated. This process is repeated with 
the offsets re-calculated at each stage until convergence to a set of 
dynamically consistent (i.e. matching in position, velocity and time) 
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perturbed conic segments is achieved. If a free-fall trajectory is 
not found, the process determines a powered trajectory with the minimum 
velocity mis -match. 

After the velocity mis-match has been minimized for the perturbed 
conic segments, the initial and final velocity offsets may be re-cal- 
culated using numerically integrated trajectory legs running between 
the same end points and times as the perturbed conic legs. At this 
point, two alternatives are possible. The first is to repeat the iter- 
ation process of the preceding steps using velocity offsets calculated 
by numerical integration at each stage.. This will provide a trajectory 
whose accuracy is limited only by the numerical precision of the inte- 
gration techniques used but will also consume a large amount of com- 
puter time. A second alternative is to accept the errors of the analytic 
technique as being well below the mid-course correction allowance and 
to use a single determination of each trajectory leg by numerical inte- 
gration as a check of the accuracy of the analytic procedure and as 
a means of determining the velocity impulse needed at each SOI entry 
and exit point to fly the trajectory predicted. 

The basic advantages of the trajectory targeting technique devel- 
oped in this, thesis are 

1) The technique is basically analytic in nature, providing 
a great reduction in the computation required. Its con- 
vergence range is wide . 

2) A continuous description of the entire trajectory is 
provided. The near -planet phases of the trajectory are 
approximated quite well by the planetocentric trajectory 
legs . 

3) By specifying the trajectory as a sequence of individual 
legs matched in position, velocity, and time,- the 


119 



determination of the trajectory is uniformly accurate 
along the trajectory. In addition, guidance objectives 
are given for each leg of the trajectory. 

4) The effects of other disturbing forces (such as non- 
gravitational effects) may be easily included in the 
perturbed conic analysis . 

5) A powered trajectory is provided for those cases which 
do not have a free-fall solution. 

The main limitation of the analytic technique lies in its accu- 
racy. The main assumption of the perturbation techniques employed is 
that each trajectory segment is basically two-body in nature. The 
presence of strong disturbing accelerations acting over extended peri- 
ods of time can cause large perturbations from the two-body reference 
legs and lead to a degradation in the accuracy of the results (as seen 
in Chapter 7). In such cases, the use of a final step employing ve- 
locity offsets calculated by numerical integration for the strongly 
perturbed legs may be necessary. 

The basic conclusion of this thesis is that the analytic targeting 
technique developed provides results sufficiently accurate for a wide 
variety of multiple swingby missions. Where its accuracy is not ade- 
quate, it may be supplemented by a final stage using numerical inte- 
gration (with the associated penalty of increased computation) to 
provide any degree of accuracy required. For heliocentric arcs (with 
the exception of those discussed in Chapter 7) , the initial and final 
velocities may be determined analytically to better than 0.4 m/sec. 
Errors m final position for a numerical integration of each leg using 
the calculated initial velocity range from 1300 km to 15,000 km with 
typical values falling in the region of 5000-8000 km. For planeto- 
centric arcs, the initial and final velocities are determined generally 
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to within 0.1 m/sec. The exceptions to this occur at Jupiter and 
Saturn for the Grand Tour example where the neglected effects of the 
disturbing accelerations due to other planets proved significant. 
Errors in final position for numerical integrations of individual legs 
using the calculated initial velocity range from 1.6 km to 240,000 km 
with the bulk of the values in the interval 2-2000 km. 

8 . 1 Contributions of the Thesis 

The author considers the following items to constitute the orig- 
inal contributions of this thesis in the field of trajectory deter- 
mination and targeting: 

1) The development of a basically analytic technique for 
the precision targeting of multiple swingby reference 
trajectories. This technique has the following features: 

aj a wide range of convergence 
b) a uniformly high level of accuracy along the 
entire trajectory 

cj specification of guidance objectives along the 
entire trajectory 

d) fast and simple to apply 

e) easily adaptable to different disturbing force 
models. 

2} Development of a method which provides an economical 
means of checking and refining the results of simple 
patched conic analyses of complicated trajectories. 

3) Development of a simple means for determining powered 
solutions for multiple swingby trajectories in cases 
where free-fall solutions do not exist. 

4} Determination of the first accurate many-body reference 
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trajectories for a multiple swingby trajectory having 
more than one intermediate swingby. 


Contributions of a secondary nature to the objective of this 
thesis are 1) development of the optimum gradient step (OGS) modifi- 
cation to the steepest descent procedure, 2) derivation of analytic 
partial derivative matrices for variations about a solution to 
Lambert's Problem, and 3) development of a perturbed conic technique 
which improves the linear range for the perturbed two -body model. 

S . 2 Recommendations for Further Study 

Several improvements and extensions of the results of this thesis 
are recommended. These are 

1) Include the effects of other disturbing accelerations 
(such as planetary oblateness, solar radiation pressure, 
low thrust, etc) in the perturbed conic model. 

2) Explore the feasibility of using the peTturbed two-body 
state transition matrix (developed in Appendix D) in the 
perturbed conic model. 

3) Study the possibility of determining (either analytically 
or numerically) the second partial derivative matrices 
for the Lambert Problem. This would allow the use of 
second order techniques to search for non-zero minima of 
the velocity mis -match cost function. 

4) Extend the perturbed conic analysis to include the deter- 
mination of injection and arrival conditions in the near- 
planet region. 

Several areas for further research using the techniques developed 
in this thesis are 
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1) Apply the analytic partial derivative matrices and so- 
lution techniques developed to the simple patched conic 
model for use in preliminary mission planning. 

2) Apply the targeting techniques developed to detailed 
mission analysis studies for the determination of launch 
windows, abort and inflight mission modification alter- 
natives, midcourse guidance requirements, and the effects 
of guidance and navigation inaccuracies. 

3) Apply the techniques for the matching of perturbed conic 
arcs to the optimization of multiple impulse orbit-to- 
orbit transfers. These transfers may be considered as 
sequences of perturbed conic coasting arcs with impulsive 
velocity changes at the matching points. The trajectory 
segment matching techniques developed here may be used 

to minimize the total impulse used for the given transfer. 
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Appendix A 

Notation and Coordinate Systems 

A, 1 Notation 

The following notation convention is used in this thesis. Examples 
are given for (3x1) vectors and (3x3) matrices but apply for any di- 
mension quantities. 

vector = column matrix 



vector transpose = row matrix 


x 


T 


[Xi x 2 x 3 ] 


vector magnitude 


x = |x| 


V 


2 2 2 
*1 + x 2 + x 3 


unit vector 


i x = X / X 

inner (or dot) product 

x • Z = x T y_ = x iyi + x 2 y 2 + x 3 y 3 
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outer (or dyadic) product 


x y 


cross product 


5. * Z “ 


x l y l 

X l y 2 

x 2 y l 

x 2 y 2 

_ X 3 y l 

X 3 y 2 


"x 2 y 3 

- x 3 y 2 

x 3^1 

’ x l?3 

*]72 

- *2^1 


x l y 3 

x 2^3 

X 3^3_ 


vector derivatives 

_ 3x. 


9x 

3a 


3a 

3x 2 

3a 

3x_ 

|_3a“ J 


3x l 

3x 1 

3x l 

3a^ 

**2 

3a 3 

^2 

3x 2 

3x 2 

3a 1 

3a 2 

3a 3 

3x 3 

3x 3 

3x 3 

3a^ 

3a 2 

^3 
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Unit Matrix 



A. 2 Coordinate Systems 

A. Z. 1 Heliocentric Coordinate System 

The heliocentric coordinate system is a non-rotating cartesian 
coordinate system centered in the sun with 

a) the positive x-axis along the line of 'inter- 
section of the earth’s equatorial plane and the 
ecliptic plane. 

b) the positive z-axis in a direction perpendicular 
to the ecliptic plane and parallel to the angular 
momentum vector of the earth about the sun 

c) the positive y-axis in the ecliptic plane located 
so as to form a right-handed coordinate system. 

A. 2 . 2 Planetocentric Coordinate System 

The planetocentric coordinate system is a non-rotating cartesian 
coordinate system centered in a planet with its axes parallel to the 
heliocentric coordinate system. 
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Appendix B 

Calculation of Conic Arcs and Their Partial Derivative Matrices 
B . 1 Calculation of Conic Trajectory Arcs 

Referring to Figure B.l, the Lambert Problem is defined as fol- 
lows : 


Given: 1) Initial position r^ and time t^ 

2) Final position and time t^ 

3) The number N of complete revolutions made about the 
central body. 

Find the initial (v^) and final (v 2 ) velocities for a two -body conic 
trajectory connecting these two points. 

The problem is solved using the following steps: 

1) Determine by some means 

a) G.^ = sgn [tt 2 - 8 2 ] (B.l) 

where 

8 = central angle traversed by the last 
incomplete revolution of the trajectory 

b) whether the trajectory is an ellipse or a 
hyperbola 

2) Calculate 

£ = r 2 - r^ , c - |£| (B.2) 

s = i (r x + r 2 + c) (B.3) 

3) For an ellipse , solve f or X in the transcendental 
equation 
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where 


Vp C t 2* t i^ = T 



(X-smX) -G^ 


(p-sxn 



and 


s(l-cos£) = (s-c) (1-cosX) 


°£ x l 2ir ; 0 < 8 _< X ’ ;0£B<.ir 

p = gravitational parameter of central body 

Calculate 

a = s/(l-cosX) 

G 2 = sgn [it 2 - X 2 ] 


For a hyperbola, solve for y in the transcendental equation 

3 

= x = ( cos g Y -i) 2 [tsinh Y -Y)-G 1 (sinh6-S)] 

where 

s(cosh<5-l) = (s-c) (coshy-1) 

and 


0< 6 < y < “ 

Calculate 

a = s/fl-coshy) 


(B.4) 
(B .5) 


(B.6) 

(B.7) 

(B.8) 

(B.9) 

(B. 10) 
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and set 


g 2 = + 1 

4) Calculate the quantities 



The initial and final velocities are given by 


^1 


v c ±C 


+ v_ 



= V i 
c — c 


V 

p -r 


[B.ll) 

(B.12) 

(B . 13) 
(B . 14) 
(B.15) 

(B.16) 

CB.17) 


where 


i r = — 
r l r l 


^ = = 2 - 
r 2 


I Ci2 


' Ii) 


(B . 18) 


The derivations for these equations may be found in [2] . 
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B.2 Calculation of Conic Partial Derivative Matrices 



This section deals with the derivation of 

matrices 







®3Ll 

3v 

1 

3v 2 


3 — 2 


3r, 

* 3r 

> 

Sr, 

9 

3r^ 


—1 


2 

—1 


—2 

and 

the partial 

derivative vectors 




3 ll 

3v 

1 

3 — 2 


^—2 


at, 

’ 3t 

9 

2 

9t l 

9 

3t 2 

The 

notation 


3a 







3x k 





3a 

3 Ik " 

V = 

3a 

3 ^k 

for 

a = 

scalar 




3a 







_ 32 kJ 




and 





“1 





3a, 

3a, 

3a, 





3x k 

3 ^k 

3z k 



3a 

7 a, “ 

3a 2 

3a 2 

3a 2 

for 


3 ^k 

v k~ 

3x k 

^k 

3z k 




3a 3 

3a 3 

3a 3 





3x k 

3 ^k 

3z k 


will 

be used for 

k = 1 

or 2 





vector 


To find the desired partial derivatives using the chain rule, it 
is necessary to derive some intermediate results. 


1 ) 


- t 


x. 2 + y v 2 + z 2 


0* 


1 , 2 


(B . 19) 


3r 


3 a 


k 

r. 


a = x, y, 
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CB.20) 


2) 


3) 


' • V k r k = — r-, 


2 _ 2 2 , „ 

" r l r 2 " 2 -1 * -2 


= r l + r 2 2 ( x l x 2 + y l y 2 + Z 1 Z 2^ 


. V 2 c = - ^ 


V = ic 


5 = 2 Cr l + r 2 + c) 


V 1 S = 2 fir 1 “ 

V 2 S = \ ( ir 2 + ic 


4) a) For an ellipse 


Then, 


3x 

3a 


(l^fix)^ l> N + ( A ' sinA ) ~ G^B-sing)] 

= |>- cos *) H - G i ( i - cose )H] 

+ J (l-cosx) 2 { (1-cosX) 2 [( 1-cosx) |-| 

J j [2 ttN + (X-smX) - G x (£-sinB) 


-s sinX 


3X 

3a. 


CB.21) 


(B.22) 


(B . 4) 


(B. 23) 
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Using 


s(l-cosg) = (s -c) (1-cos A) 


and its derivative 


33 


s sing + (1 -cos g) || = (s-c)sinA 


3A_ /_3 £ 3jc\ fl 

3a ^3a " 3a/ 1 


solve for 


(1-cosg) 



(1-cosA) 


and 


33 _ / s-c \ / sinA \ 3 A cosg-cosA 3s _ (1-cosA) 3c 

3a \ s /\sing/ 3a s sing 3a s sing 3a 


Substitute (B.4), (B.25) and (B.26) into (B.23) to get 


K-(idhx) 2 

-Gi(t£) (1-cosX) [( ¥ )(||ai)|i 


(cosg-cosA) 3s - 
s sing 3a 


— — jLS. - s sinA 
2 s 3a 1-cosA 


Groypmg terms in (B.27) 


(1 -cosA) 3c l i 
s sing 3aJj 

3a"| 

3a J 


(B.S) 

-cosA) 

(B.24) 

(B.25) 


(B.26) 


(B.27) 
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3t 

3a 



1-G ( s ~ c I iiSAl - 2 . TS i n ^ 1 3 A 

1 \ s / smgj~ 2 I-cosX| 3a 

(1-cosX) ( £ 0 . s3;cosx\ 3 t j 3s 

(l ^ os « 2 (i-ilSF)jlt O' 28 ) 


Substituting 


a 


s 

1-cosX 


CB.6) 


into (B.28) where possible 


3t 

3a 


■^(tJ isf]- 


3 at . . 
2 — sinX 


3 x [ 3s 
r s j 3a 


3X 

9a 



1 ) 3c 

sing | 3a 


(B. 29) 


Then, (B.29) may be written as 


3t 

3a 


1 3X 
Q 3a 


TT 3 S 
H 1 3a 


H 


3c 
2 3a 


(B. 30) 


where 


1 

Q 


svr[i-Gi (^) 2 


sinX ] 

smgj 


3 at , 

2 — sinX 


(B.31) 


3 35 



, T _ „ . rr- l s-c\ /cosg-cosA l 3 x 
H 1 " b l Va \~J \ s£h|i / " 2 s 


(B . 32) 


G,(s-c) 

H- = 

' Va - sing 


(B.33) 


From (B . 30) it can be seen that 


v ■ « K v - h 2 v] 


(B.34) 


since 


V^x = 0 (B.35) 

and that 

- Q C B - 36 ) 


since 


8s _ 
8x = 


0 



8x _ 
8t 


1 


b) For a hyperbola, similar procedure may be followed. 
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x 3 

« / s \ 2 

\coshy-l/ [(smhy-y) - G 1 (sinh6-5)] 


(B.8) 


Taking the derivatives, 


9t 

3a 


(coshv-l) |I - G 1 (coshfi-l) |f] 


/ s \ 2 n 

^coshy-l J L 

j^(sinhy-y) - (sinhd-S)] 


3s . , 

"3a " s sinh ^ 


lx' 

3a 


]) 


(B . 37) 


Using 


s(coshd-l) = (s-c) (coshy-1) 


(B.9) 


and its derivative 


s sinh6 |§ + (coshS-1) |f = (s-c) sinhy |J + (coshy-1) (|f - 


solve for 


(B . 38) 


(cosh<$ -1) 



(coshy-1) 


(B . 39) 


and 


36 

3a 


( s ~ c ^ /sinhy \ 3y 
\ s / \sinh6 ) 3a 


( coshy-coshS | . j)s_ 
s sinh6 ) 3a 


/ coshy-1 \ 3c 
\s sinhS / 3a 

(B. 40) 
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Substituting (B.8) , (B.39), and (B.40) into (B.37) yields 


3t 

3a 


( coshy -y {( cosh Y-l)|S - 



(coshy-1) 



+ j coshy-coshd 
I s sinhc 


> ) 3 s _ / coshy-1 J 3 cl ( 
"/3a ' Is sinhSy 3ajj 


3^ t 3s. _ s sinhy 3y 
2s 3a coshy- 1 3a 


(B.41) 


Grouping terms, 

3 

fe -{(cHSFl :)’ <^-«[ 1 - G l(¥) 2 (MS)]- 

- fa)fc°, l 'v-i) ( c ; s ^r M ) - !?}ft 

* | G i (eo3hft) t (t 2 ) <=° sI ^ - » z (5-nHw)}lf ‘ B ' 42 > 


Substituting 


a 


s 

1-coshy 


(B.10) 


into (B.42) where possible 


tiS] 

, 3 ax 
* 2 ~ 

-=inh Y ||X 


(*?) (".•SiS 0,l “) 

3 T 1 

’ 7 S J 

I 3s 
| 3a 


} 3c 



(B .43} 

(V^a” sinhd j 3a 
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Then, (B.43) may ^g wr itten as 

li-ill u 9s ,, 3c 
3a ~ Q 3a 1 3a n 2 3a 

where 



From (B.44), it may be seen th< 


V k Y ■ Q [Hj V k s - H, 7 k c] 


since 


V k t = .0 
and that 



since 


(B.44) 

(B.45) 

(B . 46) 
(B . 47) 

(B .48) 
(B.49) 
(B.50) 
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5) From (B.4) or (B.8) 


Thus 


6) a) 


Then, 


Thus, 


t = v f lT (t 2 - tj) 


(B.51) 


3a - 3a 

Tt x = 37 


3a , — 3a 

3t 2 = VTT 


For an ellipse 


a = 


1-cosX 


(B. 52) 


(B • 6) 


Vj,a = 2 [ (1-cosX) V-^s - (s sinX) V^X] 

(1-cosX) ' 


2 [(1-cosX) V^s - Qs sinX (H^V^s " 

(1-cosX) 


= (f) 2 CC I " Q sH l sin ^ V + CQsH 2 sinX;) V k c] 


V 



t p iV 


P 2 V k c] 


(B . 53) 
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where 


P, = — - OsH, sinX 
1 a 1 

P 2 = - QsH 2 sinX 

b) For a hyperbola 
s 

a 1-cosh.y 


CB.S4) 


(B.55) 


(B . 10) 


Then 


Thus 


V = 1 ; - - T7 [Cl-coshyD V k s + 


(1-coshy) 


i *- [(1-coshy) V, s + 

(l-coshy)" 2 


=(|) 2 E(| + QsH 1 sinhy) V k s 


V = (l) 2 [P lV ' P 2 V k c ^ 


s sinhy V k y] 

Qs sinhy (H 1 V k s - H 2 V k c)] 
(QsH 2 sinhy) V k c] 


(B. 56) 


where 

P 1 = I + QsH l sinhY 

P 2 = QsH 2 sinhy 


(B . 57) 


(B . 58) 
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7) For 


A 



1 

2a 


V 



1 

2A 




1 

2A 




(s-c) T 


(V v s 


V k c) 




(B.12) 


Substituting for V^a yields 


V 


■ ^ [' < V • V’ + <W - P 2 V> 


1_ 

2s' 


Collecting terms. 


^ ' ®"[( 2 J 2 ' Cs-O 2 ) v 


2s 2 (s-c) 2 


) v] 


(B . 59) 


For 


E ’ G z Vr ' it 

V - \ [ S V j - dK 
, j_ v ri - -Jj 

2B- k U 2a J 

= -ZB - [" s 2 ,V k s + ^2 V k a ] 
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= 21“ [* V k s + CP l V k s “ P 2 V k c - ) ] 
Thus, collecting terms 

V ■ - ?) V - v] 

8) For 


^ [A + B] 


(B . 13) 

VI" iv ♦ 

V k E] 


Vi j [-n- ( 

fP l 
^ 2s 2 

Cs-c) z ) * ’S" ( 2 s 2 ' 

r 1 

T 2 

1 \ , 1 ( P 2 _\lw c \ 

' [2A 

ii? 

(s-c) 2 / ZB \ 2 sVJ k j 

vf|iA 

A * i’ 


B , 

! 2A(s-c) i 2Bs 2 J k 

- r ^ 

(l 4 . 

1 \ - * 1 v, c l 

^ 4s‘ 


B / 2 ACs-c) 2 J k J 


Thus , 


Vc ’ V? [2D iV - Vt'i 


CB.61) 


where 


D 


1 



1 

4A(s-c)~^ 


1 


4Bs ‘ 


[B . 62 ) 
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D 2 = — i T 
z 45^ 


- (* * ®) ' ^ c) 


T 


Similarly, for 


Vp - VF tA - B] 

v P - V¥ tv - vi 


9) For 


Vf{ [- 


p i . k\ . 1 + J- 1 v. 

4s 2 \ A B 2ACs-c) 2 2Bs^ J 


r 


4s 


(M) 


2A(s-c) 4 J 


V 


Thus 


where 


Vp - VF I2D 3 v k s - w 1 


»3 ■ 


P 1 (1 l\ 

J? U - B ) 


1 . 1 

1 1 

4A(s-c ) c 4Bs^ 


4s 


(i-i) 


2A(s-c) ' 


l c = £ /c 


(B.63) 


(B.14) 


(B . 64) 


(B. 65) 

(B .66) 


CB.18) 
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Vic> ' 7 k 


( 4 -) 

-ji [ c, kS. • £<V> T ] 

I [ 7 k£ - ic (5 k c)T ] 


Since 


£ = L 2 - h 


2 

^2 


L Z 2 


1 

y± 


-u 


v i£ - - 1 


; V 2 c = I 


From (B.21) 


v i c = ‘ ic 


V,c = i 
2 — c 


Substituting (B.21] and (B.68) into (B.67) yields 


Vic> - - | i * lie i c T 


Vic> ’ I 1 ■ lie ie T 


Similarly, 


Vir 3 = 7 1 “ F i r lr 

i r i r 1 r 1 r ] _ r ± 


' * 2 lr 2 ^2 


(B.67) 


(B.68) 


(B.21) 


(B . 69) 


(B . 70) 
(B. 71) 
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(B . 7-2) 


Mir ) = 0 


Mir ) = 0 
^ T 1 


(B .73) 


The above intermediate 'partials may be used for the calculation 
of the desired partial derivative matrices. From 


v, <= V i + V i 
— 1 c — c P — r. 


(B.16) 


v, = V i - V n i 
—2 c — c P — r. 


(B.17) 


calculate 


V-l ■ ic t’l V c) T * Vltic) + * Vl ^! 1 


i (V-V ) T + — f- I + i i T 1 + i fV.V V 

~-c v 1 c J c L — c -cj — r^*- 1 P 


Thus 


rv V i r v i 

A’ [rj-<r] l4 ic [Vc‘<ric] 

* ±r r, v l 1 

1 L 1 p r i -r iJ 


(B i 74) 
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Similarly 


Mi ■ [I v c] 1 * ic [V c - <r ic] T ♦ i ri [»z v p] T < B - 7S > 

M 2 ■ [- I V c] 1 ♦ ic [Me * r- ic] T - ir 2 . ['i v pj T < B - 7 « 

V,V, -[-V - - vll + i [v,V - — i 1 T 

2-2 Lc c r 7 v pJ — c [2c c ^cj 


^ [Vp • £ irj 


(B . 77) 


Substitute CB.13), (B.14), CB.61), and (B.64) into (B. 74) - (B. 77) to 
get 


Ml ■ VI {[^ ce ■ B ’-| < A * E >] 1 * ic [ 2D lV - D 2 V 

* —■ J T * i, [ 2 Vl s - D 4 v lC - ^ i ] T } 


1 T 


* Vi | < E o - B i’ 1 * ic [ D i(i ri - ic> * »2ic * B i ic] T 

* ir l [Virj - V * “ 4 ic - B o irJ T J 
Ml ■ V¥[ ( E 0 - E l) 1 * < D 3 - E o) ixj irj * CW. irj *1 


+ D, i iT + (D_ + E n - D, ) i iT 1 

1 — c — r v 2 1 1 J — c — c ) 


(B.78) 
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Similarly 


where 


Define 


V 2ll " VI { E 1 1 * D 3 ir x 4 


f * 0> 3 - » 4 ) l, ij 


* ic il * 0>! - 1>2 - V ic £ } 


(B . 79) 


’1^2 = VI { ■ E 1 1 ' D 3 ir 2 irj ' C D 4 ' D 3> ir, £ 


* D 1 ic irj * (»2 * E 1 - D l> ic £ } 


(B . 80) 


^2—2 = VI {( E l - V 1 *-'&2 - D 3> lr. ir. 


2 ~ t 2 


* C » 4 - D 3 ) £ * D, ^ * ( Dl - D 2 - Ej) ^ £}' 

CB.81) 



- A + B 
1 c 



(B . 82) 

CB.83) 

CB.84) 


the matrices 
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(B.85) 



Substituting (B.85) into (B . 78) - (B. 81) yields 

V lll = VS [ CE 0 * V 1 + < D 3 - E (P M 1 + < D 4 - D 3> M 2 

+ Mg + (D 2 + E 1 - D x ) mJ (B . 86) 

^2—1 = Vf [ E 1 1 + D 3 M J + < D 3 - V M 2 + D 1 M S 

+ (D x - D 2 - M 3 ] (B.87) 

^1—2 “ Vi [■ E 1 1 ' D 3 M 4 - ( D 4 - V M 5 + D 1 M 2 

+ (D 2 + Ej. - D 1 ) M 3 J (B . 88) 

’2^2 " Vf [ (E 1 • E 2> 1 * < E 2 - “s’ H 6 * <“4 - D 3> H S 

+ D 1 Mg + CD X - D 2 - E x ) M 3 j (B.89) 

The partial derivatives with respect to time 

3 ^1 3 — 2 

3t * 3t 
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may be gotten by the chain rule using 


1) a) For an ellipse 


a 


s 

1-cosX 


(B. 6 ) 


_3a s sinX 3X 

8T " " (1-cosX) 2 3T 


a . , 3X 
= - — sinX - 5 — 
s 3 t 


Using 


3 A 
3x 


Q 


(B.36) 


3a _ _ a 2 Q 
3t s 


sinX 


CB.90) 


b) For a hyperbola 


1 -coshy 


CB.10) 


3a 

3t 


s sinhy 
( 1 -coshy ) 2 


3y 

3t 


I- sinhY u 


Using 


lx 

3t 


= Q 


(B.50) 
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2) From (B.12) - (B.1S) 


.3 A 

1 

3a 

3x 

4Aa 2 

3x 

3B 

1 

3a 

3t 

4Ba 2 

3t 

3V c 

3t 

-VII 

r 3A . 
•- — + 
L3t 

3t 

ii 

f3A 

L3t 


111 

3tJ 


Ml 

3tJ 


(B . 92) 


(B . 93) 


(B .94) 


(B.95) 


Substituting (B.90) - (B.93) 

ri ♦ ii 

3t y 2 ™ La .bJ 

Hr . jr„ ri . ii 

3t y2 La bJ 

where 

W = - sinX 
= sinhy 

• Using the above relations wi 


into (B.94) and (B.95) yields 

(B , 96) 
(B.97) 

(ellipse) 

(B . 98) 

(hyperbola) 

h (B.16) and (B.17) yields 



(B.99) 




V¥ [(i * r)ic -{l - g)ir 2 ] 


Using (B . 52) , 





3v 


2 . _ 


at. 


^2 



CB.100) 


(B.101) 


(B.102) 


B . 3 Equation Summary for Partial Derivative Matrices 


77 l “ V? [C E 0 - 

E x ) I + CD 3 - E 

0 ) M x + 

P 4 ' Dj) M 2 

* D 1 M 2* 

CD 2 - D, + E]l ) 

%] 


7 2Xl -#( E 1 1 

+ d 3 mJ - C d 4 - 

D 3> M 2 

+ D, M 5 (D 2 

* B l> M j] 




7 1—2 -V¥[- E l 1 

- D 3 M 4 - < D 4 ‘ 

V M S 

+ D 1 ^2 

* o> 2 • D i 

+ E l) M s ] 
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V 2^2 

• VF [< e i 

- e 2 ) I - 

(B 3 - E 2 ) m 6 

+ CD 4 

- V M S 


* »i “I - 

(D 3 - D x 

* *1> “ 3 ] 



M ! = 

i r ir 
r l r l 

M 3 - 

• T 

ic ic 

m 5 = 

.T 

ir 2 ic 

M 2 = 

• T 

i ri ic 

M 4 - 

. • T 

ir ir 
r 2 r l 

m 6 - 

ir £ 
r 2 r 2 

E 0 = 

A-B - 
r l 

E 1 = 

A+B 

c 

E 2 = 

A-B 

r 2 


D 


1 


D 


2 


D 


3 


D 


4 




a 

u 


+ 


+ 





1 _ l 

4A(s-c) 2 4Bs 2 

1 

2A('s-c) 2 

I + _J_ 

4A(s-c) 2 4Bs 2 

1 

2A(s-c) 2 


Ellipse 


Hyperbola 


Pi - — - QsH^ sinX 


P 1 = I + QsH l sinh Y 


P 2 = - Qs^ sinX 


?2 = QsH 2 sinhy 


H 


1 


H 


2 


r. ^nr( S-C VcosB-cosX 1 3 T 

4 va sine j ' 2 1 

G ( s ~ c ^ 1 

1 'vS" ‘ sinfS 


1 

Q 




sinX 

smB_ 


1 

Q 


o , (S-c 1 /coshy-cosh6 | 

b l v " a V~i/V sinh<S / 

G l(<7=f) siihS 

s ^[x - s,(si£) 2 g&] 


3 T 
2 s 


_3 ar 
2 s 


sinX 


+ I F sinhY 
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Figure B. 1 Lambert Problem Geometry 
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Appendix C 

Calculation of Planetary Data 
C . 1 Ephemeris Generation 

Expressions for the mean elements of the eight inner planets were 
obtained from [35] . The elements obtained from the calculations are 

£2 = longitude of ascending node 

tn = argument of perihelion 

M = mean anomaly 

L = mean longitude = £2 + to + M 

S3 = longitude of perihelion = £5 + to 

e = orbital eccentricity 

i = orbital inclination to ecliptic plane 

a = semi-major axis [in astronomical units) 

The six elements L, u, £2, e, i, a are given as expansions in the 
time in centuries 


(J.D.) - (J.D.) 

T = E. 

36525.0 

measured from the epoch Julian date for January 0.5, 1900 


(C.l) 


(J.D.) q = 2415020.0 days (C.2) 

The term (J.D.) is the current Julian date in days. The expansions, as 
determined empirically from observational astronomy are as follows . 

Mercury 

a = 0.3870984 

e = 0.20561421 + 0 . 00002046T-0 . 000000030T 2 
i = 7°0'10".37 + 6" . 699T -0".066T 2 
£2 = 47°8'45".40 + 4266". 75T + 0".626T 2 
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2 ’= 75°53 , S8".91 + 5599" . 76 T + 1".061 T 2 
L = •178°10 ' 44" . 68 + 538106654" . 80 T + 1".084 T 2 

Venus 

a = 0.72333015 

e = 0.00682069 - 0.00004774 T + 0.000000091 T 2 

V 

i = 3° 23 ' 37" . 07 + 3". 621 T - Q".Q0'3S T 2 
= 75°46'46". 73 + 3239". 46 T + 1".476 T 2 
C = 130°9'49". 8 + 5068". 93 T - 3". 515 T 2 
L = 342°46 ' 1" . 39 + 210669162". 88 T + l’M148 T 2 

Earth 

a = 1.00000013' 

e = 0.01675104 - O.OOQ0418O T - 0.000000126 T 2 
i = 0.0° 

■a = o.o 0 

a = 101°13'15".0 + 6189". 03T + 1".63 T 2 + 0".012 T 3 

L = 99 0 41'48". 04 + 129602768". 13 T + 1".089 T 2 

Mars 

a = 1.52368839 

e = 0.09331290 + 0.000092064 T - 0.000000077 T 2 

i = 1 0 51'1-". 20 - 2". 430 T + 0".0454 T 2 

= 48°47’11". 19 + 2775". 57 T - 0".005 T 2 - 0".0192 T 3 

•a = 334°13 ' 5" . 53 + 6628". 73 T + 0.4675 T 2 - 0".0043 T 3 

L = 293°44' 51". 46 f’ 68910117". 33 T + 1".1184 T 2 


Jupiter 


a = 5.202561 

e = 0.04833475 + 0.000164180 T - 0.0000004676 T 2 - 0.0000000017 T 3 
i = 1°18'31".4S - 20". 506 T + 0.014" T 2 
n = 99 0 26 ' 36" . 19' + 3637.90,8" T + 1".2680 T 2 - 0". 03064 T 3 
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ffi = 12°43 ' 15". 34 + 5795". 862 T + 3"..80258 T 2 - 
L = '238° 2 ' 57". 32 + 10930687" . 148 T + 1". 20486 

Saturn 

a = 9.554747 

e = 0.05589232 - 0.0003455 T - 0.000000728 T 2 
i = 2° 29 ' 33" .07 - 14". 108 T - 0". 05576 T 2 + 0. 
ft = 112°47 ' 25". 40 + 3143". 5025 T - 0". 54785 T 2 
(S = 91° 5 ' 53" .38 + 7050". 297 T + 2". 9749 T 2 + 0 
L = 266°33 ' 51". 76 + 4404635" . 5810 T + 1". 16835 

Uranus 

a = 19.21814 

e = 0.0463444 - 0.00002658 T + 0.000000077 T 2 
i = 0°46 ' 20" . 87 + 2". 251 T + 0".1422 T 2 
ft = 73°28 ' 37". 55 + 1795". 204 T + 4". 722 T 2 
5 = 171°32 , 55".14 + 5343". 958 T + 0”.8539 T 2 - 
L = 244°11 ' 50" . 89 + 1547508". 765 T + 1". 16835 

Neptune 

a = 30.10957 

e = 0.00899704 + 0.00000633 T - 0.000000002 T 2 
i = 1°46 ' 45" . 27 - 34". 357 T - 0".0328 T 2 
ft = 130 ° 40 ' 52" . 89 + 3956". 166 T + 0..89952 T 2 - 
S = 46°43 ' 38" . 37 + 5128". 468 T + 1". 40694 T 2 - 
L =_ 84° 27 ' 28". 78 + 791589". 291 T + 1". 15374 T 2 

The other elements are found using 


a) = a - ,0 

m = l - a 


0". 01236 T 3 
T 2 - 0.005936 T 3 


+ 0.00000000074 T 3 
00016 T 3 
- 0" . 0191 T 3 
".0166 T 3 
T 2 - 0" . 021 T 3 


0". 00218 T 3 
T 2 - 0".021 T 3 


0". 016984 T 3 
0" .002176 T 3 
- 0". 002176 T 3 


(C.3) 
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A listing of the gravitational parameter used for the planets and 
the' sun is given in Table C.l. 

Once the orbital elements of a planet have been calculated for a 
given time, its position and velocity may be calculated using the 
following relations (see [3] for derivations) . 


Ip = a 


(cos E - 1) i^+ Vap" sin E i_ 2 


Vp = Vila sin E _i^ + VpP cos E ij 


(C.4) 


where 


E = eccentric anomaly 
p = semilatus rectum 
- a(l-e 2 ) 

y = gravitational parameter of the central body (in this 
case, of the sun) 


The coordinate system defined by the unit vectors’, i^, i_ 2> I3 shown 
in Figure C.l. It can be seen that these unit vectors may be expressed 
in the ecliptic frame as 


h 


cosfl costo - sinS? sinco cos i 
sinfi cosw + cosfi sinto cos i 
sinw sin i 


—2 


cosft sinu - sinQ costa cos i 
sinO sinta + cosfi cosw cos i 
cosw sin i 


(C . 5) 
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sinfi sin i 


-3 


- cosO sin i 
cos i 


The eccentric anomaly E is related to the mean anomaly M -by Kepler's 
equation 


M = E - e sin E 


(C.6) 


This equation is transcendental and may not be solved analytically. 
The solution technique used here is a Newton-Rapheson iteration which 
sets 


M 


e sin Ei 


E k*i=V ! 


-k " ~ " 0X11 "kj 
[1 - e cos E k ] 


(C.7) 


with the initial condition 


E q = M CC.8) 

This iteration converges rapidly for orbits with small eccentricity. 

C . 2 Calculation of Sphere of Influence Radius 

In this thesis, a somewhat larger sphere of influence (SOI) is 
used than the classical Laplace SOI. This enlarged SOI is the surface 
on which the direct acceleration due to the planet equals the perturb- 
ing acceleration due to the gradient of the solar gravitational field. 
Referring to Figure C.2 it can be shown (see [3]) that an approximate 
expression (good to first order in r/£) for the acceleration on a point 
mass m due to a planet P and the sun S is 
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(C.9) 


r 


ra 


ir 


v s 


M 


3 cosa 


H 



where 


r = position of mass m with respect to the planet P 
Jl = position of the sun S with respect to the planet P 
a = angle between £■ and r 

Pg, Pp, p m = gravitational parameters of the sun, planet, 
and mass m respectively. 


Equating the direct and perturbing accelerations yields 
0.7937 < ^1 + 3 cos 2 aj 6 <1.0 


(C.10) 


Since 


and 


P m « p p or p s 


the locus of these points may be approximated by a sphere of radius 

1 

y 




Cc.iiD 


For computational purposes, l is assumed equal to the semi-major axis 
of the orbit of the planet under consideration. 

The Laplace SOI is a somewhat 'smaller surface defined approximate- 
ly fsee [3]) as 


■■ ft ) 1 


(C.12) 
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This surface is the locus of points where the ratio of 'disturbing 
acceleration to primary acceleration is equal for the equations of 
motion referred to either the sun or the planet. 

Numerical experiments by the author and by Carlson [29] show that 
a better trajectory approximation is obtained by using the enlarged 
sphere of influence. Carlson also shows that the theory of matched 
asymptotic expansion predicts an overlap of the region of validity of 
the heliocentric and planetocentric trajectory representations in the 
vicinity of this enlarged sphere of influence. A listing of the size 
of both the Laplace and the enlarged spheres of influence is given in 
Table C.2 


163 



Table C.l Gravitational Parameters 


Body 

p (km^/sec 2 ) 

Sun 

0.1327154456 (10 12 

Mercury 

0.2211924093 (10 5 ) 

Venus 

0.32528295482(10^) 

Earth 

0.39802852025(10^) 

Mars 

0.4290138858 (10 5 ) 

Jupiter 

0.12671486322(10^) 

Saturn 

0.3790137239 (10 8 ) 

Uranus 

0.580329029 (10 7 ) 

Neptune 

0.687146 34 755(10 7 ) 


Table C. 2 Sphere of Influence Radius 

Body Laplace Sphere (km) Enlarged Sphere (km) 


Mercury 

113,455.7 

318,688.4 

Venus 

616,362.0 

1,458,966.1 

Earth 

923,738.2 

2,157,378.4 

Mars 

574,520.1 

1,564,377.2 

Jupiter 

48,177,614.0 

76,638,659.8 

Saturn 

54,505,381.7 

94,129,489.7 

Uranus 

51,742,213.6 

101,287,919.7 

Neptune 

86,747,707.2 

167,883,945.9 
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ORBITAL PLANE 

periapse 
• ECLIPTIC PLANE 
ASCENDING NODE 


\lHE OF INTERSECTION 
WITH EARTH'S EQUATORIAL 


Figure C.i Coordinate System Geometry 
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Appendix D 


Properties and Application of the Two-Body State Transition Matrix 
D . 1 Properties of the Two-Body State Transition Matrix 

This section summarizes some of the more important properties of 
the two-body state transition matrix. For a more complete discussion, 
see [36] or [3] . 

The equation of motion for two -body flight is 



The variational equation for small perturbations about the nominal 
trajectory is 


Sr = G (r)6r = G [t) 6r (D.2) 

where 

Gfr) = ±3 [3 i r i r T - i] (D.3) 

r - t (t) 


Equations (D.l) and [D.2) may be put into first-order form by 
defining 


Then 


x (t) 


r (t] 
v (t) 


fix 


6r (t) 
6v (t) 



[D. 4) 
CD. 5) 


CD. 6 ) 
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and 


Sx = F(t)x 


where 


F(t) 


o r 

G (t-3 o 


The solution to (D.7) is given by 


6x (t) = $ 0 (t, t Q ) 6x(t Q ) 


CD- 7) 


CD. 8) 


CD. 9) 


where 

$g(t,t 0 ) = state transition matrix from time t Q to time t 
for the two -body trajectory. 

It may be shown [3] that the state transition matrix satisfies the 
same differential equation as the state 

« 0 Ct,t 0 ) = F(t)* 0 Ct,t 0 ) CD. 10) 

with the initial conditions 

VW = 1 CD.3.1D 

A number of analytic solutions to (D.10) in different coordinate 
systems exist (see [3], [36], [37]). The one used in this thesis is by 
Goodyear [37] , which provides a solution in generalized cartesian co- 
ordinates . 

The state transition matrix is a member of a class known as sym- 
plectic matrices. These matrices satisfy the relation 


168 



Q T JQ = J 


CD. 12) 


where 

J = 

Post-multiply (D.12) by Q”' 1 ’ and pre-raultiply by J to get 
Q' 1 = - J Q T J 

If the state transition matrix is partitioned into 


0 I 
-I 0 


Q = symplectic matrix 




then 


V Ct,t,) = 




C 0 (t,to) D Q Ct,t 0 ) 


D o T(t ’V V^'V 


■C n (t,t n ) AgCbjtQ) j 




may be calculated using (D.13). 

If (D.7) includes a disturbing term 

6x(t) = F (t) 5x(t) + f Q (t) 


where 


f Q (t) =. disturbing vector 


0 

iaCt) 


(D. 13) 


CD. 14) 


CD. IS) 


CD. 16) 
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a^(t) = disturbing acceleration evaluated along the reference 
orbit 


the solution to the perturbed motion is given by 

t 

6x(t) = * 0 Ct,t 0 ) <5x(t Q ) + j 3> 0 (t,t)f 0 (r)dx 

t 0 

See [30] for a derivation of this relation. 

D . 2 Evaluation of Perturbation Integrals by Quadrature 

In the calculation of the perturbed conic trajectory segments, 
it is necessary to evaluate the integral 


fD . 17) 


t 

s(t) * J $ 0 (t > r) f Q (t) dx CD-18) 

in (D.17) Since both terms in the integrand are known functions of time 
along the reference trajectory, the integral may be evaluated by quad- 
rature. The quadrature method chosen was Simpson’s Rule [38], which 
states that 

a+nh 

I = I u(t)dt = 5[u 0 + 4 u 2 + 2u 2 + 4 u 3 + ... 

•'a 

+ 4u _ + 2u , + 4u , + u I + R_ (D.19) 

n-3 n-2 n-1 nj n 

where 

h = step-size 

n = number of steps (even) 

Uj,= u(a+kh) 
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R n ° ■ uIV (5) 5 Ca < % < a + n h) (D.ZO) 

= truncation error 

The integral (D.18) always runs from the mid-point t^ of a tra- 
jectory segment to its end point t 2 or its initial point t ^ [see (3.12) 
and (3.13)). To take into account the fact that the integrand varys 
more rapidly with time near t^ or t^ than it does near t^, the inte- 
gral is evaluated over four sub-intervals with different step-sizes. 
These four intervals are 



This method insures that the interval for 1^ is divided into 
an even number of steps approximately d in length. 


For a small planet (Mercury, Venus, Earth, Mars) the parameter T g is 
assigned the values 


T s = 0.5 days (heliocentric leg) 

T g = 0.02 days (planetocentric leg) (D.21) 
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while for a large planet (Jupiter, Saturn, Uranus, Neptune) it has 
the' values 


T = 2,5 days (heliocentric leg) 

(D.22) 

T g = 0.5 days (planetocentric leg) 

These values are picked to provide a balance between running time and 
accuracy . 

The value of the integral in (D.18) is given by the sum of the 
integrals over the sub - intervals . 

1 = I 1 + I 2 + T 3 + J 4 (D. 23) 

To illustrate the accuracy of this procedure, a comparison between 
the perturbations obtained by quadrature and those found by numerical 
integration is given for 

1) Heliocentric Ellipse (Earth to Venus ) 

t ± -= 2441478.8 days t M = 2441556.4599 days 

. t 2 = 2441634.11977 days 

Quadrature Numerical Integration Difference 

SrC^) 35,164.67 km 35,267.35 km 102.78 km 

dv(t^) 39.723 m/sec 40.237 m/sec .514 m/sec 

2) Planetocentric Hyperbola (about Venus) 

t x = 2441634.11977 days t M = 2441636.0597 days 

t 2 - 2441637.99955 days 
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luadrature 


Numerical Integration Difference 


9 


<$r(t^) 1482.11 km 1512.57 km 

Svft^ 26.470 m/sec 26.489 m/sec 


3^ Heliocentric Ellipse (Earth to Jupiter) 


30.46 km 
.190 m/sec 


t x = 2443787.0 days t M = 2444040.0151 days 

= 2444293.03027 days 

Quadrature Numerical Integration Difference 



Errors in the calculation of the perturbations by the analytic 
quadrature method are on the order of 1% or less. Position perturba- 
tions tend to be more accurate than the velocity perturbations. This 
is due to the fact that the true A and B sub -matrices in the state 
transition matrix (see (d.14)) are more accurately approximated by 
their value on the two-body reference trajectory than are the C and D 
sub-matrices . 

D . 3 Calculation of Perturbed State Transition Matrix 

It was stated in (D.7) and (D.10) that the state variation Sx and 
the state transition matrix both satisfy the same differential equation. 
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5x(t) = F(t)6x(t) 


CD. 7) 


*Ct,-Jt 0 ) = F(t)*(t,t n ) CD. 10) 

i U 

Furthermore, if a perturbing term f is added to (D.7), its solution 
becomes 


for 


<5x(t) 


$Ct,t Q )6x(t) 


+ 


j 4> (t ,x)f (t) dr 
r 0 


CD. 17) 


6x(t) = FCt)5x(t) + f(t) (D. 16) 

A similar form may be derived for the perturbed state transition 
matrix. Let 


* 0 (t»t 0 ) = F 0 (t) * 0 Ct,t 0 ) {D. 24) 

be the differential equation for the pure two-body problem while 

*Ct,t 0 ) = F(t)O(t,t 0 ) CD.2S) 

is the differential equation for the perturbed two -body problem. Note 
that 

F(t) = F 0 (t) + F d (t) 



1 

l-t 

0 

1 


o 

0 

1 



+ 


. 

_G 0. 


- G d 


(D. 26) 
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where 


9a, 
G d = S*- 


(D. 27) 


Substituting (D.26) into (D.25) yields 


i’Ct.tg) = Fg (t) 0 (t , tg) + F d (t)$ Ct,tg) 


(D.28) 


4>(t,tg) = §g(t,tg) + 6$Ct,t 0 ) 


CD . 293 


Then,. (D.28) becomes 


(D . 29) 


* 0 (t,t 0 ) + ««>(t,t n ) = F n (t)* ft (t,t n ) + F n (t)6$(t,t n ) 


Q J ~ • * Q < 


+ F d (t)$ 0 (t,t 0 ) + F d (t)6$(t,t 0 ) (D. 30) 


Eliminating the two-body terms and neglecting the product F d (t) 6$ (t , tg) 
leaves 


<S*.(t,t 0 ) = FgC^d^-Ct.tf,) + F^ (t) $ n (t ,t n ) 


O' T ^d vt -' ,v 0 v ’ 0 ' 


(D.31) 


By comparing (D.31) with (D.16), the solution analogous to (D.17) is 


5$(t,t 0 ) = $ 0 (t,t 0 )6$(t 0 . 


L. 

.t 0 ) + f $ o (t,T ^ F d^' r ' ) dT 

J t n (D.32) 


Inserting (D.32) into (D.29) and recognizing that 


5§ C t n ,t f) ) = 0 


yields 
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« , 0 (t,x)F d (T)$ 0 ^ T » t o) <i ' i: 


(D , 33) 


= $ ft (t, 


V * / 


which is an approximate solution for a many-body state transition 
matrix. 


The errors associated with the use of the pure two -body state 
transition matrix rather than the actual many body matrix were the 
major source of inaccuracy in the perturbed conic analysis. However, 
since these errors were not excessive, it was felt that the increased 
computation time associated with the above calculation of an approxi- 
mate many-body state transition matrix was not justified. No numerical 
studies on the. accuracy of the above technique were performed. 
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Appendix E 


Trajectory Description Data 
E. 1 Description of Tables 

This, appendix presents the detailed data needed for the specifi- 
cation of -the trajectories determined in Chapter 5-7. The following 
sections contain tables of 

1) Planetary position in the heliocentric ecliptic coordinate 
system at each sphere of influence (SOI) entry and exit 
time . 

2) Planetary velocity in the heliocentric ecliptic coordinate 
system at each SOI entry and exit time. 

3) Position of SOI entry and exit points in the planetocen- 
tric ecliptic coordinate frame. 

4) True heliocentric velocity at each SOI entry and exit 
point . 

5) True planetocentric velocity at each SOI entry and exit 
point . 

The true velocities at each SOI entry and exit point are deter- 
mined by numerical integration of each trajectory leg as described in 
Section 5.4. The time of passage through each SOI entry and exit point 
may be found in Chapters 5-7. Planets are referred to by number in the 
tables in the following manner 


Number 

Planet 

Number 

Planet 

1 

Mercury 

5 

Jupiter 

2 

Venus 

6 

Saturn 

3 

Earth 

7 

Uranus 

4 

Mars 

8 

Neptune 
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E . 2 Dual Planet Reconnaissance Trajectory 


Planetary Position fkm) 




Point 

Planet 

X 

Z 

z_ 

1 

3 

-27,683,569 

-149,351,633 

0 

2 

2 

-85,946,304 

64,291,179 

5,849,706 

3 

2 

-92,494,950 

54,464,844 

6,090,584 

4 

4 

9,298,184 

-216,541,101 

-4,783,633 

5 

4 

20,235,530 

-214,917,173 

-5,017,472 

6 

3 

1-50,071,922 

1,341,038 

0 

Planetary Velocity Qcm/sec) 



Point 

Planet 

V 

X 

Iz 

^z 

1 

3 

28.803850 

-5.536705 

o 

o 

2 

2 

-21.115118 • 

-28.219406 

0.824872 

3 

2 

-17.916849 

-30.347939 

0.610865 

4 

4 

25.142150 

3.112230 

-0.550617 

5 

4 

25.057758 

4.344065 

-0.522634 

6 

3 

-0.752478 

29.680227 

O 

O 

Planetocentric 

Coordinates of SOI Entry 

and Exit Points 

flan) 

Point 

Planet 

X 

z 

z 

1 

3 

-1,618,847 

1,228,772 

-723,696 

2 

2 

1,289,792 

-678,612 

-67,123 

3 

2 

-1,457,233 

56,086 

43,682 

4 

4 

1,139,936 

1,065,458 

112,352 

5 

4 

-1,151,504 

-1,053,602 

1,060,012 

6 

3 

2,124,226 

309,891 

214,270 
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True Heliocentric Velocity at SOI Entry and Exit .Points ucm/secj 


‘ Point Planet 




!x 


v 


z 


1 

2 

3 

4 

5 

6 


3 
2 
2 

4 
4 
3 


25.411804 

-28.702396 

-26.543009 

19.925562 

19.783087 

-11.562989 


-3.308020 

-24.100878 

-30.113072 

- 1.771693 

-0.479968 

27.887921 


-1.496488 

1.225602 

0.865247 

- 1.114716 

0.012257 

0.623309 


True Planetocentric Velocity at SOI Entry and Exit Points (km/sec) 


Point Planet 




v 


z 


1 

2 

3 

4 

5 

6 


3 

- 

- 

- 

2 

-7.587318 

4.118546 , 

0.400726 

2 

-8.626212 

. 0.234868 

0.254384 

4 

-5.216591 

-4.883929 

-0.564099 

4 

-5.274674 

-4.824036 

0.534891 

3 

. 

.. 

- 


E.3 Grand Tour Trajectory 
Planetary Position [km] 

Point Planet x X. 


1 

2 

3 

4 

5 

6 

7 

8 


3 

146,373,093 

30,840,195 

0 

5 

-730,966,280 

341,268,726 

14,992,492' 

5 

-793,050,422 

172,374,356 

17,066,856 

6 

-1,392,750,100 

-367,498,888 

61,896,089 

6 

-1,351,245,940 

-523,607,941 

62,947,475 

7 

-508,598,411 

-2,816,199,860 

-3,939,139 

.7 

-417,872,638 

-2,835,394,500 

-5,185,687 

8 

825,783,251 

-4,452,305,730 

72,635,448 
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Planetary Velocity (km/sec) 


Point 

Planet 

V 

X 

Li 

V 

z 

1 

3 

-6.627671 

29.040716 

0.0 

2 

S 

-5.683434 

-11.237136 

0.172861 

3 

5 

-2.929249 

-12.166462 

0.114945 

4 

6 

1.929444 

-9.356882 

0.08S318 

5 

6 

2.952499 

-9.024942 

0.038810 

6 

7 

6.653486 

-1.521178 

-0.091945 

7 

7 

6.689117 

-1.304052 

-0.091596 

8 

8 

5.299626 

1.020036 

-0.143199 

Planetocentric Coordinates of SOI Entry 

and Exit Points 

(km) 

Point 

Planet 

X 

Z 

2_ 

1 

3 

-422,358 

2,090,193 

327,087 

2 

5 

52,179,576 

-56,112,733 

1,461,922 

3 

S 

-76,427,535 

-88,025 

5,684,043 

4 

6 

90,303,036 

25,960,699 

5,635,961 

5 

6 

18,853,800 

-92,130,386 

4,10.9,294 

6 

7 

-7,355,051 

100,977,635 

2,943,333 

7. 

7 

36,980,520 

-94,149,082 

5,256,825 

8 

8 

-68,924,481 

152,861,419 

8,235,400 

True Heliocentric Velocity at SOI Entry 

and Exit Points 

(km/, sec) 

Point 

Planet 


Zz 


1 

3 

-8.729043 

38.773437 

1.658246 

2’ 

5 

: 12. 546516 

-3.447447 

0 . 060104 

3 

5 

-13.280658 

-12.488954 

0.910924 

4' 

.6 

-8.652141 

-12.357576 

0.746294 

5 

6 

5.181211 

-19.749870 

-0.442043 

6 

7 

7.716297 

-16.446444 

-0.533078 
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Point 


Planet 


7 7 12.174013 -15.208984 0.691257 

8 8 12.136149 -14.246205 0.684301 


True Planetocentric Velocity at SOI Entry and Exit Points (km/sec) 


Point 

Planet 

v x 

Iz 

V 

z 

1 

3 

- 

- 

- 

2 

5 

-6.862569 

7.789369 

-0.232962 

3 

5 

-10.351326 

-0.322564 

0.795978 

4 

6 

-10.581829 

-3.000683 

0.660993 

5 

6 

2.228758 

-10.724730 

-0.480851 

6 

7 

1.062830 

-14.925182 

-0.441132 

7 

7 

5.484978 

-13.904982 

0.782847 

8 

8 

- 

- 

- 

E.4 Periodic 

Trajectory 



Planetary Position (km) 




Point 

Planet 

X 

y 

z 

1 

3 

115,517,131 

-98,137,483 

0 

2 

2 

-98,406,259 

42,926,751 

6,271,195 

3 

2 

104,690,651'. 

24,078,614 

6,371,004 

4 

2 

-97,993,845 

43,854,474 

6,260,237 

5 

2 

104,439,345 

25,127,936 

6,371,140 

6 

2 

-97,465,377 

45,010,278 

6,245,750 

7 

2 

104,152,073 

26,274,532 

6,370,544 

8 

3 

149,582,514 

10,210,861 

0 

9 

3 

144,544,696 

38,009,310 

0 

10 

3 

149,721,363 

8, 463, -429 

0 

11 

3 

145,020,497 

36,278,351 

0 

12 

2 

-96,870,027 

-47,526,925 

4,924,934 
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Planetary Velocity (km/sec) 


Point 

Planet 

^x 

Zx 

V 

Z 

1 

3 

18.800505 

22.595294 

0.0 

2 

2 

-14.161646 

-32.263762 

0.367108 

3 

2 

-8.029999 

-34.286904 

-0.014566 

4 

2 

-14.463523 

-32.130379 

0.386568 

5 

2 

-8.371178 

-34.206640 

0.006423 

6 

2 

-14.839637 

-31.959383 

0.410832 

7 

2 

-8.744011 

-34.114777 

0.029399 

8 

3 

-2.515094 

29.612395 

0.0 

9 

3 

-8.062106 

28.702144 

0.0 

10 

3 

-2.167525 

29.633935 

0.0 

11 

3 

-7.715572 

28.790893 

O 

o 

12 

2 

15.182219 

-31.602193 

-1.314630 

Planetocentric 

Coordinates of SOI 

Entry and Exit Points fkm) 

Point 

Planet 

X 

X 

z 

1 

3 

247,388 

-639,293 

-2,045,577 

2 

2 

726,703 

511,959 

-1,156,885 

3 

2 

-1,279,966 

527,150 

460,850 

4 

2 

1,283,436 

r S47,542 

-426,113 

5 

2 

-770,819 

372,484 

-1,181,388 

6 

2 

748,191 

-351,377 

1,202,217 

7 

2 

9,682 

-767,274 

-1,240,878 

8 

3 

337,614 

1,396,334 

-1,609,518 

9 

3 

1,614,003 

107,287 

1,427,503 

10 

3 

1,637,219 

-113,046 

-1,400,362 

11 

3 

1,107,649 

-1,140,039 

-1,458,666 

12 

2 

t 538 ,663 

296,754 

-1,323,012 
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True Heliocentric Velocity at SOI Entry 

and Exit Points (km/sec) 

Point 

Planet 


V 

JL 


1 

3 

17.323915 

20.157965 

-3.267526 

2 

2 

-16.616974 

-34.071864 

4.356735 

3 

2 

-12.484461 

-32.404149 

1.506887 

4 

2 

-18.920016 

-30.231828 

1.908059 

5 

2 

-11.046420 

-32.921844 

-4.047927 

6 

2 

-17.512140 

-30.658423 

-3.641063 

7 

2 

-8.650968 

-36.809919 

-4.202766 

8 

3 

-3.262805 

26.747620 

3.2S4240 

9 

3 

-4.699129 

28.972881 

2.858207 

10 

3 

1.194127 

29.890394 

2.857028 

11 

3 

-5.477269 

26.475428 

-2.992208 

12 

2 

16.261390 

-33.864871 

3.969588 

True Planetocentric 

Velocity at SOI Entry and Exit 

Points (km/sec' 

Point 

Planet 



V 

Z 

1 

3 

- 

- 

- 

2 

2 

-2.455124 

-1.807893' 

3.989942 

3 

2 

-4.455058 

1.883022 

1.531741 

4 

2 

-4.457173 

1.899298 

1.532338 

5 

2 

-2.631615 

1.269900 

-4.102211 

6 

2 

-2.625721 

1.284945 

-4.097597 

7 

2 

0.093025 

-2.694363 

-4.232504 

8 

3 

-0.746902 

-2.862650 

3.253536 

9 

3 

3.351399 

0.267886 

2.882569 

10 

3 

3.348385 

0.253327 

2.879679 

11 

3 

2.239589 

-2.317191 

-2.991828 

12 

2 

_ 

_ 

_ 
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